Level set: composition (II)

include/adaflo/level_set_okz_advance_concentration.h

@@ -101,15 +101,13 @@ class LevelSetOKZSolverAdvanceConcentration
     const VectorType &                            vel_solution,
     const VectorType &                            vel_solution_old,
     const VectorType &                            vel_solution_old_old,
-    const double &                                global_omega_diameter,
     const AlignedVector<VectorizedArray<double>> &cell_diameters,
     const AffineConstraints<double> &             constraints,
     const ConditionalOStream &                    pcout,
     const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor<dim> &boundary,
     const MatrixFree<dim> &                                             matrix_free,
     const LevelSetOKZSolverAdvanceConcentrationParameter &              parameters,
-    double &                              global_max_velocity,
-    const DiagonalPreconditioner<double> &preconditioner);
+    const DiagonalPreconditioner<double> &                              preconditioner);
 
   virtual void
   advance_concentration(const double dt);
@@ -167,11 +165,11 @@ class LevelSetOKZSolverAdvanceConcentration
   /**
    * Physics section
    */
-  const double &                                global_omega_diameter;           // [i]
+  double                                        global_omega_diameter;           // [i]
   const AlignedVector<VectorizedArray<double>> &cell_diameters;                  // [i]
   const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor<dim> boundary;   // [i]
   AlignedVector<VectorizedArray<double>>                 artificial_viscosities; // [-]
-  double &                                               global_max_velocity;    // [o]
+  double                                                 global_max_velocity;    // [o]
   AlignedVector<Tensor<1, dim, VectorizedArray<double>>> evaluated_convection;   // [o]
 
   /**

include/adaflo/util.h

@@ -18,6 +18,8 @@
 
 #include <deal.II/distributed/tria.h>
 
+#include <deal.II/matrix_free/matrix_free.h>
+
 using namespace dealii;
 
 /**
@@ -52,4 +54,78 @@ locally_owned_subdomain(const MeshType &mesh)
   return tria_parallel != nullptr ? tria_parallel->locally_owned_subdomain() : 0;
 }
 
+template <int dim>
+void
+compute_cell_diameters(const MatrixFree<dim, double> &         matrix_free,
+                       const unsigned int                      dof_index,
+                       AlignedVector<VectorizedArray<double>> &cell_diameters,
+                       double &                                cell_diameter_min,
+                       double &                                cell_diameter_max)
+{
+  cell_diameters.resize(matrix_free.n_cell_batches());
+
+  cell_diameter_min = std::numeric_limits<double>::max();
+  cell_diameter_max = 0.0;
+
+  // to find the cell diameters, we compute the maximum and minimum eigenvalue
+  // of the Jacobian transformation from the unit to the real cell. We check
+  // all face centers and the center of the cell and take the respective
+  // minimum and maximum there to cover most of the cell geometry
+  std::vector<Point<dim>> face_centers;
+  {
+    Point<dim> center;
+    for (unsigned int d = 0; d < dim; ++d)
+      center[d] = 0.5;
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        Point<dim> p1 = center;
+        p1[d]         = 0;
+        face_centers.push_back(p1);
+        p1[d] = 1.;
+        face_centers.push_back(p1);
+      }
+    face_centers.push_back(center);
+  }
+
+  const auto &dof_handler   = matrix_free.get_dof_handler(dof_index);
+  const auto &triangulation = dof_handler.get_triangulation();
+
+  LAPACKFullMatrix<double> mat(dim, dim);
+  FEValues<dim>            fe_values(*matrix_free.get_mapping_info().mapping,
+                          dof_handler.get_fe(),
+                          Quadrature<dim>(face_centers),
+                          update_jacobians);
+  for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell)
+    {
+      VectorizedArray<double> diameter = VectorizedArray<double>();
+      for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(cell); ++v)
+        {
+          typename DoFHandler<dim>::active_cell_iterator dcell =
+            matrix_free.get_cell_iterator(cell, v, 1);
+          fe_values.reinit(dcell);
+          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
+            {
+              mat = 0;
+              for (unsigned int d = 0; d < dim; ++d)
+                for (unsigned int e = 0; e < dim; ++e)
+                  mat(d, e) = fe_values.jacobian(q)[d][e];
+              mat.compute_eigenvalues();
+              for (unsigned int d = 0; d < dim; ++d)
+                {
+                  diameter[v] = std::max(diameter[v], std::abs(mat.eigenvalue(d)));
+                  cell_diameter_min =
+                    std::min(cell_diameter_min, std::abs(mat.eigenvalue(d)));
+                }
+            }
+          if (1U + dcell->level() == triangulation.n_global_levels())
+            cell_diameter_max = std::max(diameter[v], cell_diameter_max);
+        }
+      cell_diameters[cell] = diameter;
+    }
+  cell_diameter_min =
+    -Utilities::MPI::max(-cell_diameter_min, get_communicator(triangulation));
+  cell_diameter_max =
+    Utilities::MPI::max(cell_diameter_max, get_communicator(triangulation));
+}
+
 #endif

source/level_set_okz.cc

@@ -188,14 +188,12 @@ LevelSetOKZSolver<dim>::LevelSetOKZSolver(const FlowParameters &parameters_in,
         this->navier_stokes.solution.block(0),
         this->navier_stokes.solution_old.block(0),
         this->navier_stokes.solution_old_old.block(0),
-        this->global_omega_diameter,
         this->cell_diameters,
         this->constraints,
         this->pcout,
         bcs,
         this->matrix_free,
         params,
-        this->global_max_velocity,
         this->preconditioner);
   }
 }

source/level_set_okz_advance_concentration.cc

@@ -64,6 +64,36 @@ namespace
 
     return Utilities::MPI::max(max_velocity, get_communicator(dof_handler));
   }
+
+  template <int dim, int spacedim>
+  double
+  diameter_on_coarse_grid(const Triangulation<dim, spacedim> &tria)
+  {
+    const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
+    std::vector<bool>                   boundary_vertices(vertices.size(), false);
+
+    for (const auto &cell : tria.cell_iterators_on_level(0))
+      for (const unsigned int face : cell->face_indices())
+        if (cell->face(face)->at_boundary())
+          for (unsigned int i = 0; i < cell->face(face)->n_vertices(); ++i)
+            boundary_vertices[cell->face(face)->vertex_index(i)] = true;
+
+    // now traverse the list of boundary vertices and check distances.
+    // since distances are symmetric, we only have to check one half
+    double                            max_distance_sqr = 0;
+    std::vector<bool>::const_iterator pi               = boundary_vertices.begin();
+    const unsigned int                N                = boundary_vertices.size();
+    for (unsigned int i = 0; i < N; ++i, ++pi)
+      {
+        std::vector<bool>::const_iterator pj = pi + 1;
+        for (unsigned int j = i + 1; j < N; ++j, ++pj)
+          if ((*pi == true) && (*pj == true) &&
+              ((vertices[i] - vertices[j]).norm_square() > max_distance_sqr))
+            max_distance_sqr = (vertices[i] - vertices[j]).norm_square();
+      }
+
+    return std::sqrt(max_distance_sqr);
+  }
 } // namespace
 
 
@@ -154,14 +184,12 @@ LevelSetOKZSolverAdvanceConcentration<dim>::LevelSetOKZSolverAdvanceConcentratio
   const VectorType &                            vel_solution,
   const VectorType &                            vel_solution_old,
   const VectorType &                            vel_solution_old_old,
-  const double &                                global_omega_diameter,
   const AlignedVector<VectorizedArray<double>> &cell_diameters,
   const AffineConstraints<double> &             constraints,
   const ConditionalOStream &                    pcout,
   const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor<dim> &boundary,
   const MatrixFree<dim> &                                             matrix_free,
   const LevelSetOKZSolverAdvanceConcentrationParameter &              parameters,
-  double &                                                            global_max_velocity,
   const DiagonalPreconditioner<double> &                              preconditioner)
   : parameters(parameters)
   , solution(solution)
@@ -176,10 +204,9 @@ LevelSetOKZSolverAdvanceConcentration<dim>::LevelSetOKZSolverAdvanceConcentratio
   , constraints(constraints)
   , pcout(pcout)
   , time_stepping(parameters.time)
-  , global_omega_diameter(global_omega_diameter)
+  , global_omega_diameter(0.0)
   , cell_diameters(cell_diameters)
   , boundary(boundary)
-  , global_max_velocity(global_max_velocity)
   , preconditioner(preconditioner)
 {}
 
@@ -474,6 +501,10 @@ LevelSetOKZSolverAdvanceConcentration<dim>::advance_concentration(const double d
   this->time_stepping.set_time_step(dt);
   this->time_stepping.next();
 
+  if (global_omega_diameter == 0.0)
+    global_omega_diameter =
+      diameter_on_coarse_grid(matrix_free.get_dof_handler().get_triangulation());
+
   if (evaluated_convection.size() !=
       this->matrix_free.n_cell_batches() *
         this->matrix_free.get_n_q_points(parameters.quad_index))

source/level_set_okz_preconditioner.cc

@@ -178,6 +178,17 @@ initialize_projection_matrix(
     ilu_projection_matrix.initialize(projection_matrix);
   }
 }
+template void
+initialize_projection_matrix<1, double, VectorizedArray<double>>(
+  const MatrixFree<1, double, VectorizedArray<double>> &matrix_free,
+  const AffineConstraints<double> &                     constraints_normals,
+  const unsigned int                                    dof_index,
+  const unsigned int                                    quad_index,
+  const double &                                        epsilon_used,
+  const double &                                        epsilon,
+  const AlignedVector<VectorizedArray<double>> &        cell_diameters,
+  BlockMatrixExtension &                                projection_matrix,
+  BlockILUExtension &                                   ilu_projection_matrix);
 
 template void
 initialize_projection_matrix<2, double, VectorizedArray<double>>(

source/two_phase_base.cc

@@ -279,66 +279,11 @@ TwoPhaseBaseAlgorithm<dim>::initialize_data_structures()
 
   // find relevant epsilon for smoothing by taking the largest mesh size of
   // cells close to the interface (here: cells on finest level)
-  epsilon_used        = 0;
-  minimal_edge_length = global_omega_diameter;
-  cell_diameters.resize(this->matrix_free.n_cell_batches());
-
-  // to find the cell diameters, we compute the maximum and minimum eigenvalue
-  // of the Jacobian transformation from the unit to the real cell. We check
-  // all face centers and the center of the cell and take the respective
-  // minimum and maximum there to cover most of the cell geometry
-  std::vector<Point<dim>> face_centers;
-  {
-    Point<dim> center;
-    for (unsigned int d = 0; d < dim; ++d)
-      center[d] = 0.5;
-    for (unsigned int d = 0; d < dim; ++d)
-      {
-        Point<dim> p1 = center;
-        p1[d]         = 0;
-        face_centers.push_back(p1);
-        p1[d] = 1.;
-        face_centers.push_back(p1);
-      }
-    face_centers.push_back(center);
-  }
-  LAPACKFullMatrix<double> mat(dim, dim);
-  FEValues<dim>            fe_values(this->mapping,
-                          navier_stokes.get_fe_p(),
-                          Quadrature<dim>(face_centers),
-                          update_jacobians);
-  for (unsigned int cell = 0; cell < this->matrix_free.n_cell_batches(); ++cell)
-    {
-      VectorizedArray<double> diameter = VectorizedArray<double>();
-      for (unsigned int v = 0;
-           v < this->matrix_free.n_active_entries_per_cell_batch(cell);
-           ++v)
-        {
-          typename DoFHandler<dim>::active_cell_iterator dcell =
-            this->matrix_free.get_cell_iterator(cell, v, 1);
-          fe_values.reinit(dcell);
-          for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
-            {
-              mat = 0;
-              for (unsigned int d = 0; d < dim; ++d)
-                for (unsigned int e = 0; e < dim; ++e)
-                  mat(d, e) = fe_values.jacobian(q)[d][e];
-              mat.compute_eigenvalues();
-              for (unsigned int d = 0; d < dim; ++d)
-                {
-                  diameter[v] = std::max(diameter[v], std::abs(mat.eigenvalue(d)));
-                  minimal_edge_length =
-                    std::min(minimal_edge_length, std::abs(mat.eigenvalue(d)));
-                }
-            }
-          if (1U + dcell->level() == this->triangulation.n_global_levels())
-            epsilon_used = std::max(diameter[v], epsilon_used);
-        }
-      cell_diameters[cell] = diameter;
-    }
-  minimal_edge_length =
-    -Utilities::MPI::max(-minimal_edge_length, get_communicator(triangulation));
-  epsilon_used = Utilities::MPI::max(epsilon_used, get_communicator(triangulation));
+  compute_cell_diameters(this->matrix_free,
+                         1 /*dof_index_p*/,
+                         cell_diameters,
+                         minimal_edge_length,
+                         epsilon_used);
 
   this->pcout << "Mesh size (largest/smallest element length at finest level): "
               << epsilon_used << " / " << minimal_edge_length << std::endl;