Merge pull request #26 from peterrum/ls_composition_2

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/level_set_okz_compute_curvature.h

@@ -23,7 +23,6 @@
 
 #include <adaflo/block_matrix_extension.h>
 #include <adaflo/diagonal_preconditioner.h>
-#include <adaflo/level_set_okz_compute_normal.h>
 #include <adaflo/navier_stokes.h>
 #include <adaflo/time_stepping.h>
 #include <adaflo/util.h>
@@ -76,7 +75,6 @@ class LevelSetOKZSolverComputeCurvature
 {
 public:
   LevelSetOKZSolverComputeCurvature(
-    LevelSetOKZSolverComputeNormal<dim> &                  normal_operator,
     const AlignedVector<VectorizedArray<double>> &         cell_diameters,
     const LinearAlgebra::distributed::BlockVector<double> &normal_vector_field,
     const AffineConstraints<double> &                      constraints_curvature,
@@ -120,11 +118,6 @@ class LevelSetOKZSolverComputeCurvature
    */
   const LevelSetOKZSolverComputeCurvatureParameter parameters; // [i]
 
-  /**
-   * Other operators.
-   */
-  LevelSetOKZSolverComputeNormal<dim> &normal_operator; // [i]
-
   /**
    * Vector section
    */

include/adaflo/level_set_okz_compute_normal.h

@@ -56,8 +56,11 @@ struct LevelSetOKZSolverComputeNormalParameter
   /**
    * TODO
    */
-
   bool approximate_projections;
+  /**
+   * TODO
+   */
+  double damping_scale_factor = 4.0;
 };
 
 template <int dim>
@@ -70,7 +73,7 @@ class LevelSetOKZSolverComputeNormal
   LevelSetOKZSolverComputeNormal(
     BlockVectorType &                              normal_vector_field,
     BlockVectorType &                              normal_vector_rhs,
-    VectorType &                                   solution,
+    const VectorType &                             level_set_field,
     const AlignedVector<VectorizedArray<double>> & cell_diameters,
     const double &                                 epsilon_used,
     const double &                                 minimal_edge_length,
@@ -112,7 +115,7 @@ class LevelSetOKZSolverComputeNormal
    */
   BlockVectorType & normal_vector_field; // [o]
   BlockVectorType & normal_vector_rhs;   // [-]
-  const VectorType &vel_solution;        // [i]
+  const VectorType &level_set_solution;  // [i]
 
   /**
    * MatrixFree

include/adaflo/level_set_okz_reinitialization.h

@@ -68,7 +68,6 @@ class LevelSetOKZSolverReinitialization
   using BlockVectorType = LinearAlgebra::distributed::BlockVector<double>;
 
   LevelSetOKZSolverReinitialization(
-    LevelSetOKZSolverComputeNormal<dim> &             normal_operator,
     const BlockVectorType &                           normal_vector_field,
     const AlignedVector<VectorizedArray<double>> &    cell_diameters,
     const double &                                    epsilon_used,
@@ -84,7 +83,6 @@ class LevelSetOKZSolverReinitialization
     bool &                                            first_reinit_step,
     const MatrixFree<dim, double> &                   matrix_free)
     : parameters(parameters)
-    , normal_operator(normal_operator)
     , solution(solution)
     , solution_update(solution_update)
     , system_rhs(system_rhs)
@@ -102,11 +100,11 @@ class LevelSetOKZSolverReinitialization
   {}
 
   // performs reinitialization
-  virtual void
-  reinitialize(const double       dt,
-               const unsigned int stab_steps,
-               const unsigned int diff_steps                              = 0,
-               const bool         diffuse_cells_with_large_curvature_only = false);
+  void
+  reinitialize(const double                     dt,
+               const unsigned int               stab_steps,
+               const unsigned int               diff_steps     = 0,
+               const std::function<void(bool)> &compute_normal = [](const bool) {});
 
   void
   reinitialization_vmult(VectorType &      dst,
@@ -133,11 +131,6 @@ class LevelSetOKZSolverReinitialization
    */
   const LevelSetOKZSolverReinitializationParameter parameters;
 
-  /**
-   * Other operators.
-   */
-  LevelSetOKZSolverComputeNormal<dim> &normal_operator;
-
   /**
    * Vector section
    */

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, dof_index);
+          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

@@ -105,7 +105,6 @@ LevelSetOKZSolver<dim>::LevelSetOKZSolver(const FlowParameters &parameters_in,
     params.time.time_step_size_min   = this->parameters.time_step_size_min;
 
     this->reinit_operator = std::make_unique<LevelSetOKZSolverReinitialization<dim>>(
-      *this->normal_operator,
       this->normal_vector_field,
       this->cell_diameters,
       this->epsilon_used,
@@ -133,7 +132,6 @@ LevelSetOKZSolver<dim>::LevelSetOKZSolver(const FlowParameters &parameters_in,
     params.curvature_correction    = this->parameters.curvature_correction;
 
     this->curvatur_operator = std::make_unique<LevelSetOKZSolverComputeCurvature<dim>>(
-      *this->normal_operator,
       this->cell_diameters,
       this->normal_vector_field,
       this->constraints_curvature,
@@ -188,14 +186,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);
   }
 }
@@ -464,8 +460,16 @@ template <int dim>
 void
 LevelSetOKZSolver<dim>::compute_curvature(const bool diffuse_large_values)
 {
-  TimerOutput::Scope timer(*this->timer, "LS compute curvature.");
-  curvatur_operator->compute_curvature(diffuse_large_values);
+  // This function computes the curvature from the normal field. Could also
+  // compute the curvature directly from C, but that is less accurate. TODO:
+  // include that variant by a parameter
+  {
+    this->compute_normal(false);
+  }
+  {
+    TimerOutput::Scope timer(*this->timer, "LS compute curvature.");
+    curvatur_operator->compute_curvature(diffuse_large_values);
+  }
 }
 
 
@@ -543,11 +547,15 @@ LevelSetOKZSolver<dim>::reinitialize(const unsigned int stab_steps,
                                      const unsigned int diff_steps,
                                      const bool diffuse_cells_with_large_curvature_only)
 {
+  (void)diffuse_cells_with_large_curvature_only;
+
   TimerOutput::Scope timer(*this->timer, "LS reinitialization step.");
   reinit_operator->reinitialize(this->time_stepping.step_size(),
                                 stab_steps,
                                 diff_steps,
-                                diffuse_cells_with_large_curvature_only);
+                                [&](const bool fast_computation) {
+                                  this->compute_normal(fast_computation);
+                                });
 }
 
 

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))