Add CH output, in progress

chaos-polymtl · Mar 26, 2024 · 955ac1e · 955ac1e
1 parent d80099d
commit 955ac1e
Show file tree

Hide file tree

Showing 2 changed files with 80 additions and 8 deletions.
diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h
@@ -315,14 +315,7 @@ class CahnHilliard : public AuxiliaryPhysics<dim, GlobalVectorType>
    * @param[in] display_precision Number of outputted digits.
    */
   void
-  output_newton_update_norms(const unsigned int display_precision) override
-  {
-    this->pcout << std::setprecision(display_precision)
-                << "\t||dx||_L2 = " << std::setw(6) << newton_update.l2_norm()
-                << std::setw(6)
-                << "\t||dx||_Linfty = " << std::setprecision(display_precision)
-                << newton_update.linfty_norm() << std::endl;
-  }
+  output_newton_update_norms(const unsigned int display_precision) override;
 
 
 private:

diff --git a/source/solvers/cahn_hilliard.cc b/source/solvers/cahn_hilliard.cc
@@ -1322,6 +1322,85 @@ CahnHilliard<dim>::apply_phase_filter()
     }
 }
 
+template <int dim>
+void
+CahnHilliard<dim>::output_newton_update_norms(
+  const unsigned int display_precision)
+{
+      auto mpi_communicator = triangulation->get_communicator();
+
+      FEValuesExtractors::Scalar phase_order(0);
+      FEValuesExtractors::Scalar chemical_potential(1);
+
+      ComponentMask phase_order_mask = fe->component_mask(phase_order);
+      ComponentMask chemical_potential_mask = fe->component_mask(chemical_potential);
+
+      const std::vector<IndexSet> index_set_phase_order =
+        DoFTools::locally_owned_dofs_per_component(dof_handler, phase_order_mask);
+      const std::vector<IndexSet> index_set_chemical_potential =
+        DoFTools::locally_owned_dofs_per_component(dof_handler, chemical_potential_mask);
+
+      double local_sum = 0.0;
+      double local_max = std::numeric_limits<double>::lowest();
+
+
+          for (auto j = index_set_phase_order[0].begin();
+               j != index_set_phase_order[0].end();
+               j++)
+            {
+              double dof_newton_update = newton_update[*j];
+
+              local_sum += dof_newton_update * dof_newton_update;
+
+              local_max = std::max(local_max, dof_newton_update);
+            }
+
+
+      double global_phase_order_l2_norm =
+        std::sqrt(Utilities::MPI::sum(local_sum, mpi_communicator));
+      double global_phase_order_linfty_norm =
+        Utilities::MPI::max(local_max, mpi_communicator);
+
+      local_sum = 0.0;
+      local_max = std::numeric_limits<double>::lowest();
+
+      for (auto j = index_set_chemical_potential[1].begin();
+           j != index_set_chemical_potential[1].end();
+           j++)
+        {
+          double dof_newton_update = newton_update[*j];
+
+          local_sum += dof_newton_update * dof_newton_update;
+
+          local_max = std::max(local_max, dof_newton_update);
+        }
+
+      double global_chemical_potential_l2_norm =
+        std::sqrt(Utilities::MPI::sum(local_sum, mpi_communicator));
+      double global_chemical_potential_linfty_norm =
+        Utilities::MPI::max(local_max, mpi_communicator);
+
+      this->pcout << std::setprecision(display_precision)
+                  << "\n\t||du||_L2 = " << std::setw(6)
+                  << newton_update.l2_norm() << std::setw(6)
+                  << "\t||du||_Linfty = "
+                  << std::setprecision(display_precision)
+                  << newton_update.linfty_norm() << std::endl;
+
+      this->pcout << std::setprecision(display_precision)
+                  << "\n\t||dphi||_L2 = " << std::setw(6)
+                  << global_phase_order_l2_norm << std::setw(6)
+                  << "\t||dphi||_Linfty = "
+                  << std::setprecision(display_precision)
+                  << global_phase_order_linfty_norm << std::endl;
+      this->pcout << std::setprecision(display_precision)
+                  << "\t||deta||_L2 = " << std::setw(6) << global_chemical_potential_l2_norm
+                  << std::setw(6) << "\t||deta||_Linfty = "
+                  << std::setprecision(display_precision)
+                  << global_chemical_potential_linfty_norm << std::endl;
+
+}
+
 template std::pair<Tensor<1, 2>, Tensor<1, 2>>
 CahnHilliard<2>::calculate_barycenter<GlobalVectorType>(
   const GlobalVectorType &solution,