From b79b3c5b4f530bc7a3f56c8d32462dcff7a90d0c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 15 Dec 2020 06:50:06 +0100 Subject: [PATCH 01/16] Make LevelSetOKZSolverAdvanceConcentration::global_max_velocity local --- include/adaflo/level_set_okz_advance_concentration.h | 5 ++--- source/level_set_okz.cc | 1 - source/level_set_okz_advance_concentration.cc | 2 -- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/include/adaflo/level_set_okz_advance_concentration.h b/include/adaflo/level_set_okz_advance_concentration.h index 886da539..50dcf7bd 100644 --- a/include/adaflo/level_set_okz_advance_concentration.h +++ b/include/adaflo/level_set_okz_advance_concentration.h @@ -108,8 +108,7 @@ class LevelSetOKZSolverAdvanceConcentration const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor &boundary, const MatrixFree & matrix_free, const LevelSetOKZSolverAdvanceConcentrationParameter & parameters, - double & global_max_velocity, - const DiagonalPreconditioner &preconditioner); + const DiagonalPreconditioner & preconditioner); virtual void advance_concentration(const double dt); @@ -171,7 +170,7 @@ class LevelSetOKZSolverAdvanceConcentration const AlignedVector> &cell_diameters; // [i] const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor boundary; // [i] AlignedVector> artificial_viscosities; // [-] - double & global_max_velocity; // [o] + double global_max_velocity; // [o] AlignedVector>> evaluated_convection; // [o] /** diff --git a/source/level_set_okz.cc b/source/level_set_okz.cc index 8ac96a48..17690835 100644 --- a/source/level_set_okz.cc +++ b/source/level_set_okz.cc @@ -195,7 +195,6 @@ LevelSetOKZSolver::LevelSetOKZSolver(const FlowParameters ¶meters_in, bcs, this->matrix_free, params, - this->global_max_velocity, this->preconditioner); } } diff --git a/source/level_set_okz_advance_concentration.cc b/source/level_set_okz_advance_concentration.cc index 4b5a7045..3a32ce66 100644 --- a/source/level_set_okz_advance_concentration.cc +++ b/source/level_set_okz_advance_concentration.cc @@ -161,7 +161,6 @@ LevelSetOKZSolverAdvanceConcentration::LevelSetOKZSolverAdvanceConcentratio const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor &boundary, const MatrixFree & matrix_free, const LevelSetOKZSolverAdvanceConcentrationParameter & parameters, - double & global_max_velocity, const DiagonalPreconditioner & preconditioner) : parameters(parameters) , solution(solution) @@ -179,7 +178,6 @@ LevelSetOKZSolverAdvanceConcentration::LevelSetOKZSolverAdvanceConcentratio , global_omega_diameter(global_omega_diameter) , cell_diameters(cell_diameters) , boundary(boundary) - , global_max_velocity(global_max_velocity) , preconditioner(preconditioner) {} From 240da41f5b339c1f34f725657236e8bb19ccd0a6 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 15 Dec 2020 07:48:47 +0100 Subject: [PATCH 02/16] Compute global_omega_diameter within LevelSetOKZSolverAdvanceConcentration --- .../level_set_okz_advance_concentration.h | 3 +- source/level_set_okz.cc | 1 - source/level_set_okz_advance_concentration.cc | 37 ++++++++++++++++++- 3 files changed, 36 insertions(+), 5 deletions(-) diff --git a/include/adaflo/level_set_okz_advance_concentration.h b/include/adaflo/level_set_okz_advance_concentration.h index 50dcf7bd..63985ad3 100644 --- a/include/adaflo/level_set_okz_advance_concentration.h +++ b/include/adaflo/level_set_okz_advance_concentration.h @@ -101,7 +101,6 @@ class LevelSetOKZSolverAdvanceConcentration const VectorType & vel_solution, const VectorType & vel_solution_old, const VectorType & vel_solution_old_old, - const double & global_omega_diameter, const AlignedVector> &cell_diameters, const AffineConstraints & constraints, const ConditionalOStream & pcout, @@ -166,7 +165,7 @@ class LevelSetOKZSolverAdvanceConcentration /** * Physics section */ - const double & global_omega_diameter; // [i] + double global_omega_diameter; // [i] const AlignedVector> &cell_diameters; // [i] const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor boundary; // [i] AlignedVector> artificial_viscosities; // [-] diff --git a/source/level_set_okz.cc b/source/level_set_okz.cc index 17690835..1674be75 100644 --- a/source/level_set_okz.cc +++ b/source/level_set_okz.cc @@ -188,7 +188,6 @@ LevelSetOKZSolver::LevelSetOKZSolver(const FlowParameters ¶meters_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, diff --git a/source/level_set_okz_advance_concentration.cc b/source/level_set_okz_advance_concentration.cc index 3a32ce66..c981fdc0 100644 --- a/source/level_set_okz_advance_concentration.cc +++ b/source/level_set_okz_advance_concentration.cc @@ -64,6 +64,36 @@ namespace return Utilities::MPI::max(max_velocity, get_communicator(dof_handler)); } + + template + double + diameter_on_coarse_grid(const Triangulation &tria) + { + const std::vector> &vertices = tria.get_vertices(); + std::vector 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::const_iterator pi = boundary_vertices.begin(); + const unsigned int N = boundary_vertices.size(); + for (unsigned int i = 0; i < N; ++i, ++pi) + { + std::vector::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,7 +184,6 @@ LevelSetOKZSolverAdvanceConcentration::LevelSetOKZSolverAdvanceConcentratio const VectorType & vel_solution, const VectorType & vel_solution_old, const VectorType & vel_solution_old_old, - const double & global_omega_diameter, const AlignedVector> &cell_diameters, const AffineConstraints & constraints, const ConditionalOStream & pcout, @@ -175,7 +204,7 @@ LevelSetOKZSolverAdvanceConcentration::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) , preconditioner(preconditioner) @@ -472,6 +501,10 @@ LevelSetOKZSolverAdvanceConcentration::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)) From e5dc601b2399fa78e14655acd5a8acbef8cd1240 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 15 Dec 2020 10:08:44 +0100 Subject: [PATCH 03/16] Create utility function compute_cell_diameters() --- include/adaflo/util.h | 76 ++++++++++++++++++++++++++++++++++++++++ source/two_phase_base.cc | 65 +++------------------------------- 2 files changed, 81 insertions(+), 60 deletions(-) diff --git a/include/adaflo/util.h b/include/adaflo/util.h index f1e276b3..246037bf 100644 --- a/include/adaflo/util.h +++ b/include/adaflo/util.h @@ -18,6 +18,8 @@ #include +#include + using namespace dealii; /** @@ -52,4 +54,78 @@ locally_owned_subdomain(const MeshType &mesh) return tria_parallel != nullptr ? tria_parallel->locally_owned_subdomain() : 0; } +template +void +compute_cell_diameters(const MatrixFree & matrix_free, + const unsigned int dof_index, + AlignedVector> &cell_diameters, + double & cell_diameter_min, + double & cell_diameter_max) +{ + cell_diameters.resize(matrix_free.n_cell_batches()); + + cell_diameter_min = std::numeric_limits::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> face_centers; + { + Point center; + for (unsigned int d = 0; d < dim; ++d) + center[d] = 0.5; + for (unsigned int d = 0; d < dim; ++d) + { + Point 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 mat(dim, dim); + FEValues fe_values(*matrix_free.get_mapping_info().mapping, + dof_handler.get_fe(), + Quadrature(face_centers), + update_jacobians); + for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell) + { + VectorizedArray diameter = VectorizedArray(); + for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(cell); ++v) + { + typename DoFHandler::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 diff --git a/source/two_phase_base.cc b/source/two_phase_base.cc index 13f65b76..22aa3fef 100644 --- a/source/two_phase_base.cc +++ b/source/two_phase_base.cc @@ -279,66 +279,11 @@ TwoPhaseBaseAlgorithm::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> face_centers; - { - Point center; - for (unsigned int d = 0; d < dim; ++d) - center[d] = 0.5; - for (unsigned int d = 0; d < dim; ++d) - { - Point p1 = center; - p1[d] = 0; - face_centers.push_back(p1); - p1[d] = 1.; - face_centers.push_back(p1); - } - face_centers.push_back(center); - } - LAPACKFullMatrix mat(dim, dim); - FEValues fe_values(this->mapping, - navier_stokes.get_fe_p(), - Quadrature(face_centers), - update_jacobians); - for (unsigned int cell = 0; cell < this->matrix_free.n_cell_batches(); ++cell) - { - VectorizedArray diameter = VectorizedArray(); - for (unsigned int v = 0; - v < this->matrix_free.n_active_entries_per_cell_batch(cell); - ++v) - { - typename DoFHandler::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; From 7496f51c3f77c5fc59d668435073b9ca1307cc0c Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 11:00:40 +0100 Subject: [PATCH 04/16] 1d projection added --- source/level_set_okz_preconditioner.cc | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/source/level_set_okz_preconditioner.cc b/source/level_set_okz_preconditioner.cc index cfabb8aa..80264dae 100644 --- a/source/level_set_okz_preconditioner.cc +++ b/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>( + const MatrixFree<1, double, VectorizedArray> &matrix_free, + const AffineConstraints & constraints_normals, + const unsigned int dof_index, + const unsigned int quad_index, + const double & epsilon_used, + const double & epsilon, + const AlignedVector> & cell_diameters, + BlockMatrixExtension & projection_matrix, + BlockILUExtension & ilu_projection_matrix); template void initialize_projection_matrix<2, double, VectorizedArray>( From 5228896bb406c5d211028168a3b0582decd994f9 Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 11:57:26 +0100 Subject: [PATCH 05/16] bug fix in compute_cell_diameters --- include/adaflo/util.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/adaflo/util.h b/include/adaflo/util.h index 246037bf..56b3d06f 100644 --- a/include/adaflo/util.h +++ b/include/adaflo/util.h @@ -101,7 +101,7 @@ compute_cell_diameters(const MatrixFree & matrix_free, for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(cell); ++v) { typename DoFHandler::active_cell_iterator dcell = - matrix_free.get_cell_iterator(cell, v, 1); + 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) { From 44dc62da5da04db0425a913d8de1cdcf64822c15 Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 11:57:55 +0100 Subject: [PATCH 06/16] compute normal: remove hardcoded indices --- source/level_set_okz_compute_normal.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/level_set_okz_compute_normal.cc b/source/level_set_okz_compute_normal.cc index 4612cf65..59c5e9e2 100644 --- a/source/level_set_okz_compute_normal.cc +++ b/source/level_set_okz_compute_normal.cc @@ -89,7 +89,7 @@ LevelSetOKZSolverComputeNormal::local_compute_normal( // The second input argument below refers to which constrains should be used, // 4 means constraints_normals const unsigned int n_q_points = ls_degree == -1 ? 0 : 2 * ls_degree; - FEEvaluation phi(data, 4, 2); + FEEvaluation phi(data, parameters.dof_index_normal, parameters.quad_index); const VectorizedArray min_diameter = make_vectorized_array(this->epsilon_used / this->parameters.epsilon); // cast avoids compile errors, but we always use the path without casting @@ -128,8 +128,8 @@ LevelSetOKZSolverComputeNormal::local_compute_normal_rhs( // The second input argument below refers to which constrains should be used, // 4 means constraints_normals and 2 means constraints (for LS-function) const unsigned int n_q_points = ls_degree == -1 ? 0 : 2 * ls_degree; - FEEvaluation normal_values(data, 4, 2); - FEEvaluation ls_values(data, 2, 2); + FEEvaluation normal_values(data, parameters.dof_index_normal, parameters.quad_index); + FEEvaluation ls_values(data, parameters.dof_index_ls, parameters.quad_index); for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { From f3cb70ae884a7a0bbf96d4901ed83a3fd236a700 Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 12:27:46 +0100 Subject: [PATCH 07/16] ComputeNormal: parameter damping_scale_factor added --- include/adaflo/level_set_okz_compute_normal.h | 5 ++++- source/level_set_okz_compute_normal.cc | 14 ++++++++++---- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/include/adaflo/level_set_okz_compute_normal.h b/include/adaflo/level_set_okz_compute_normal.h index 003a39ad..39456e8a 100644 --- a/include/adaflo/level_set_okz_compute_normal.h +++ b/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 diff --git a/source/level_set_okz_compute_normal.cc b/source/level_set_okz_compute_normal.cc index 59c5e9e2..ece2a684 100644 --- a/source/level_set_okz_compute_normal.cc +++ b/source/level_set_okz_compute_normal.cc @@ -89,7 +89,9 @@ LevelSetOKZSolverComputeNormal::local_compute_normal( // The second input argument below refers to which constrains should be used, // 4 means constraints_normals const unsigned int n_q_points = ls_degree == -1 ? 0 : 2 * ls_degree; - FEEvaluation phi(data, parameters.dof_index_normal, parameters.quad_index); + FEEvaluation phi(data, + parameters.dof_index_normal, + parameters.quad_index); const VectorizedArray min_diameter = make_vectorized_array(this->epsilon_used / this->parameters.epsilon); // cast avoids compile errors, but we always use the path without casting @@ -101,7 +103,7 @@ LevelSetOKZSolverComputeNormal::local_compute_normal( phi.read_dof_values(src); phi.evaluate(true, true); const VectorizedArray damping = - Number(4.) * + Number(parameters.damping_scale_factor) * Utilities::fixed_power<2>( std::max(min_diameter, cell_diameters[cell] / static_cast(ls_degree))); for (unsigned int q = 0; q < phi.n_q_points; ++q) @@ -128,8 +130,12 @@ LevelSetOKZSolverComputeNormal::local_compute_normal_rhs( // The second input argument below refers to which constrains should be used, // 4 means constraints_normals and 2 means constraints (for LS-function) const unsigned int n_q_points = ls_degree == -1 ? 0 : 2 * ls_degree; - FEEvaluation normal_values(data, parameters.dof_index_normal, parameters.quad_index); - FEEvaluation ls_values(data, parameters.dof_index_ls, parameters.quad_index); + FEEvaluation normal_values(data, + parameters.dof_index_normal, + parameters.quad_index); + FEEvaluation ls_values(data, + parameters.dof_index_ls, + parameters.quad_index); for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { From e996eb2c1c33c4b1f654018e9ac1f7608e7d2a5e Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 14:33:02 +0100 Subject: [PATCH 08/16] 1d instance of solution curvature added --- source/level_set_okz_compute_curvature.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/source/level_set_okz_compute_curvature.cc b/source/level_set_okz_compute_curvature.cc index 12b8c92b..7182c040 100644 --- a/source/level_set_okz_compute_curvature.cc +++ b/source/level_set_okz_compute_curvature.cc @@ -332,5 +332,6 @@ LevelSetOKZSolverComputeCurvature::compute_curvature(const bool) +template class LevelSetOKZSolverComputeCurvature<1>; template class LevelSetOKZSolverComputeCurvature<2>; template class LevelSetOKZSolverComputeCurvature<3>; From 97ee84bdef33c7d600a3be2ba59c33e76afe45f8 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 15 Dec 2020 15:45:07 +0100 Subject: [PATCH 09/16] Add specialization for 1D --- source/level_set_okz_compute_curvature.cc | 93 +++++++++++++++++++---- 1 file changed, 79 insertions(+), 14 deletions(-) diff --git a/source/level_set_okz_compute_curvature.cc b/source/level_set_okz_compute_curvature.cc index 7182c040..46f068ff 100644 --- a/source/level_set_okz_compute_curvature.cc +++ b/source/level_set_okz_compute_curvature.cc @@ -124,6 +124,81 @@ LevelSetOKZSolverComputeCurvature::local_compute_curvature( } } +namespace +{ + template + VectorizedArray + normalize(const VectorizedArray &in, bool &all_zero) + { + VectorizedArray vec; + + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + vec[v] = in[v] >= 1e-2 ? 1.0 : (in[v] <= -1e-2 ? -1.0 : 0.0); + + all_zero = false; // TODO? + + return vec; + } + + template + static Tensor<1, dim, VectorizedArray> + normalize(const Tensor<1, dim, VectorizedArray> &normal, bool &all_zero) + { + Tensor<1, dim, VectorizedArray> vec = normal; + + const VectorizedArray normal_norm = normal.norm(); + for (unsigned int d = 0; d < VectorizedArray::size(); ++d) + if (normal_norm[d] > 1e-2) + { + all_zero = false; + for (unsigned int e = 0; e < dim; ++e) + vec[e][d] /= normal_norm[d]; + } + else + for (unsigned int e = 0; e < dim; ++e) + vec[e][d] = 0; + + return vec; + } + + template + static VectorizedArray + get_divergence(const FEEvaluation &phi, + const unsigned int q) + { + return phi.get_divergence(q); + } + + template + static VectorizedArray + get_divergence(const FEEvaluation<1, + fe_degree, + n_q_points_1d, + n_components_, + Number, + VectorizedArrayType> &phi, + const unsigned int q) + { + (void)q; + (void)phi; + return 0.0; + } +} // namespace + template @@ -159,26 +234,16 @@ LevelSetOKZSolverComputeCurvature::local_compute_curvature_rhs( bool all_zero = true; for (unsigned int i = 0; i < normal_values.dofs_per_component; ++i) { - Tensor<1, dim, VectorizedArray> normal = normal_values.get_dof_value(i); - const VectorizedArray normal_norm = normal.norm(); - for (unsigned int d = 0; d < VectorizedArray::size(); ++d) - if (normal_norm[d] > 1e-2) - { - all_zero = false; - for (unsigned int e = 0; e < dim; ++e) - normal[e][d] /= normal_norm[d]; - } - else - for (unsigned int e = 0; e < dim; ++e) - normal[e][d] = 0; - normal_values.submit_dof_value(normal, i); + normal_values.submit_dof_value(normalize(normal_values.get_dof_value(i), + all_zero), + i); } if (all_zero == false) { normal_values.evaluate(false, true, false); for (unsigned int q = 0; q < normal_values.n_q_points; ++q) - curv_values.submit_value(-normal_values.get_divergence(q), q); + curv_values.submit_value(-get_divergence(normal_values, q), q); curv_values.integrate(true, false); curv_values.distribute_local_to_global(dst); } From 3438aff7c2cc8b0085ceafb597bb20209cdafed4 Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 15:47:44 +0100 Subject: [PATCH 10/16] make variable names more verbose --- source/level_set_okz_compute_normal.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/level_set_okz_compute_normal.cc b/source/level_set_okz_compute_normal.cc index ece2a684..4016b626 100644 --- a/source/level_set_okz_compute_normal.cc +++ b/source/level_set_okz_compute_normal.cc @@ -52,7 +52,7 @@ template LevelSetOKZSolverComputeNormal::LevelSetOKZSolverComputeNormal( BlockVectorType & normal_vector_field, BlockVectorType & normal_vector_rhs, - VectorType & solution, + const VectorType & level_set_field, const AlignedVector> & cell_diameters, const double & epsilon_used, const double & minimal_edge_length, @@ -65,7 +65,7 @@ LevelSetOKZSolverComputeNormal::LevelSetOKZSolverComputeNormal( : parameters(parameters) , normal_vector_field(normal_vector_field) , normal_vector_rhs(normal_vector_rhs) - , vel_solution(solution) + , level_set_solution(level_set_field) , matrix_free(matrix_free) , constraints_normals(constraints_normals) , cell_diameters(cell_diameters) @@ -142,7 +142,7 @@ LevelSetOKZSolverComputeNormal::local_compute_normal_rhs( normal_values.reinit(cell); ls_values.reinit(cell); - ls_values.read_dof_values_plain(this->vel_solution); + ls_values.read_dof_values_plain(this->level_set_solution); ls_values.evaluate(false, true, false); for (unsigned int q = 0; q < normal_values.n_q_points; ++q) @@ -218,7 +218,7 @@ LevelSetOKZSolverComputeNormal::compute_normal(const bool fast_computation) &LevelSetOKZSolverComputeNormal::template local_compute_normal_rhs, \ this, \ this->normal_vector_rhs, \ - this->vel_solution) + this->level_set_solution) EXPAND_OPERATIONS(OPERATION); #undef OPERATION From 91fe4a44a26e41fadced3b4991f42752814bbd39 Mon Sep 17 00:00:00 2001 From: Magdalena Karoline Schreter Date: Tue, 15 Dec 2020 15:49:14 +0100 Subject: [PATCH 11/16] make variable names more verbose --- include/adaflo/level_set_okz_compute_normal.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/adaflo/level_set_okz_compute_normal.h b/include/adaflo/level_set_okz_compute_normal.h index 39456e8a..e6769b53 100644 --- a/include/adaflo/level_set_okz_compute_normal.h +++ b/include/adaflo/level_set_okz_compute_normal.h @@ -73,7 +73,7 @@ class LevelSetOKZSolverComputeNormal LevelSetOKZSolverComputeNormal( BlockVectorType & normal_vector_field, BlockVectorType & normal_vector_rhs, - VectorType & solution, + const VectorType & level_set_field, const AlignedVector> & cell_diameters, const double & epsilon_used, const double & minimal_edge_length, @@ -115,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 From 91abb2046220620a641f5ee9675fc4138b418b44 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 15 Dec 2020 16:00:04 +0100 Subject: [PATCH 12/16] Instantiate LevelSetOKZSolverReinitialization for 1D --- source/level_set_okz_reinitialization.cc | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/source/level_set_okz_reinitialization.cc b/source/level_set_okz_reinitialization.cc index 77650cf5..5fbd07a4 100644 --- a/source/level_set_okz_reinitialization.cc +++ b/source/level_set_okz_reinitialization.cc @@ -105,6 +105,24 @@ LevelSetOKZSolverReinitialization::local_reinitialize( } +namespace +{ + template + VectorizedArray + normalize(const VectorizedArray &in) + { + return std::abs(in); + } + + template + static VectorizedArray + normalize(const Tensor<1, dim, VectorizedArray> &in) + { + return in.norm(); + } +} // namespace + + template template @@ -147,8 +165,8 @@ LevelSetOKZSolverReinitialization::local_reinitialize_rhs( Tensor<1, dim, VectorizedArray> grad = phi.get_gradient(q); if (first_reinit_step) { - Tensor<1, dim, VectorizedArray> normal = normals.get_value(q); - normal /= std::max(make_vectorized_array(1e-4), normal.norm()); + auto normal = normals.get_value(q); + normal /= std::max(make_vectorized_array(1e-4), normalize(normal)); evaluated_normal[cell * phi.n_q_points + q] = normal; } // take normal as it was for the first reinit step @@ -346,5 +364,6 @@ LevelSetOKZSolverReinitialization::reinitialize(const double dt, this->time_stepping.next(); } +template class LevelSetOKZSolverReinitialization<1>; template class LevelSetOKZSolverReinitialization<2>; template class LevelSetOKZSolverReinitialization<3>; From 62cf3543b91b2430b7857b4055e334d3c6e22ae0 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 16 Dec 2020 08:06:40 +0100 Subject: [PATCH 13/16] Eliminate normal_operator from LevelSetOKZSolverComputeCurvature --- include/adaflo/level_set_okz_compute_curvature.h | 7 ------- source/level_set_okz.cc | 14 +++++++++++--- source/level_set_okz_compute_curvature.cc | 7 ------- 3 files changed, 11 insertions(+), 17 deletions(-) diff --git a/include/adaflo/level_set_okz_compute_curvature.h b/include/adaflo/level_set_okz_compute_curvature.h index 9916b36f..08d1350e 100644 --- a/include/adaflo/level_set_okz_compute_curvature.h +++ b/include/adaflo/level_set_okz_compute_curvature.h @@ -23,7 +23,6 @@ #include #include -#include #include #include #include @@ -76,7 +75,6 @@ class LevelSetOKZSolverComputeCurvature { public: LevelSetOKZSolverComputeCurvature( - LevelSetOKZSolverComputeNormal & normal_operator, const AlignedVector> & cell_diameters, const LinearAlgebra::distributed::BlockVector &normal_vector_field, const AffineConstraints & constraints_curvature, @@ -120,11 +118,6 @@ class LevelSetOKZSolverComputeCurvature */ const LevelSetOKZSolverComputeCurvatureParameter parameters; // [i] - /** - * Other operators. - */ - LevelSetOKZSolverComputeNormal &normal_operator; // [i] - /** * Vector section */ diff --git a/source/level_set_okz.cc b/source/level_set_okz.cc index 1674be75..9cc4ff97 100644 --- a/source/level_set_okz.cc +++ b/source/level_set_okz.cc @@ -133,7 +133,6 @@ LevelSetOKZSolver::LevelSetOKZSolver(const FlowParameters ¶meters_in, params.curvature_correction = this->parameters.curvature_correction; this->curvatur_operator = std::make_unique>( - *this->normal_operator, this->cell_diameters, this->normal_vector_field, this->constraints_curvature, @@ -462,8 +461,17 @@ template void LevelSetOKZSolver::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 + { + TimerOutput::Scope timer(*this->timer, "LS compute normal."); + this->compute_normal(false); + } + { + TimerOutput::Scope timer(*this->timer, "LS compute curvature."); + curvatur_operator->compute_curvature(diffuse_large_values); + } } diff --git a/source/level_set_okz_compute_curvature.cc b/source/level_set_okz_compute_curvature.cc index 46f068ff..4f790cda 100644 --- a/source/level_set_okz_compute_curvature.cc +++ b/source/level_set_okz_compute_curvature.cc @@ -50,7 +50,6 @@ template LevelSetOKZSolverComputeCurvature::LevelSetOKZSolverComputeCurvature( - LevelSetOKZSolverComputeNormal & normal_operator, const AlignedVector> & cell_diameters, const LinearAlgebra::distributed::BlockVector &normal_vector_field, const AffineConstraints & constraints_curvature, @@ -65,7 +64,6 @@ LevelSetOKZSolverComputeCurvature::LevelSetOKZSolverComputeCurvature( std::shared_ptr & projection_matrix, std::shared_ptr & ilu_projection_matrix) : parameters(parameters) - , normal_operator(normal_operator) , solution_curvature(solution_curvature) , rhs(system_rhs) , solution_ls(solution_ls) @@ -327,11 +325,6 @@ template void LevelSetOKZSolverComputeCurvature::compute_curvature(const bool) { - // 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 - normal_operator.compute_normal(false); - // compute right hand side rhs = 0; From 6ee1fb5c9698212fe3d5bf9dd0286657a819e37d Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 16 Dec 2020 08:55:21 +0100 Subject: [PATCH 14/16] Eliminate normal_operator from LevelSetOKZSolverReinitialization --- include/adaflo/level_set_okz_reinitialization.h | 15 ++++----------- source/level_set_okz.cc | 7 +++++-- source/level_set_okz_reinitialization.cc | 11 ++++++----- 3 files changed, 15 insertions(+), 18 deletions(-) diff --git a/include/adaflo/level_set_okz_reinitialization.h b/include/adaflo/level_set_okz_reinitialization.h index 686af095..a862e839 100644 --- a/include/adaflo/level_set_okz_reinitialization.h +++ b/include/adaflo/level_set_okz_reinitialization.h @@ -68,7 +68,6 @@ class LevelSetOKZSolverReinitialization using BlockVectorType = LinearAlgebra::distributed::BlockVector; LevelSetOKZSolverReinitialization( - LevelSetOKZSolverComputeNormal & normal_operator, const BlockVectorType & normal_vector_field, const AlignedVector> & cell_diameters, const double & epsilon_used, @@ -84,7 +83,6 @@ class LevelSetOKZSolverReinitialization bool & first_reinit_step, const MatrixFree & matrix_free) : parameters(parameters) - , normal_operator(normal_operator) , solution(solution) , solution_update(solution_update) , system_rhs(system_rhs) @@ -103,10 +101,10 @@ 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); + reinitialize(const double dt, + const unsigned int stab_steps, + const unsigned int diff_steps = 0, + const std::function &compute_normal = [](const bool) {}); void reinitialization_vmult(VectorType & dst, @@ -133,11 +131,6 @@ class LevelSetOKZSolverReinitialization */ const LevelSetOKZSolverReinitializationParameter parameters; - /** - * Other operators. - */ - LevelSetOKZSolverComputeNormal &normal_operator; - /** * Vector section */ diff --git a/source/level_set_okz.cc b/source/level_set_okz.cc index 9cc4ff97..3ed6cba1 100644 --- a/source/level_set_okz.cc +++ b/source/level_set_okz.cc @@ -105,7 +105,6 @@ LevelSetOKZSolver::LevelSetOKZSolver(const FlowParameters ¶meters_in, params.time.time_step_size_min = this->parameters.time_step_size_min; this->reinit_operator = std::make_unique>( - *this->normal_operator, this->normal_vector_field, this->cell_diameters, this->epsilon_used, @@ -549,11 +548,15 @@ LevelSetOKZSolver::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); + }); } diff --git a/source/level_set_okz_reinitialization.cc b/source/level_set_okz_reinitialization.cc index 5fbd07a4..18061bab 100644 --- a/source/level_set_okz_reinitialization.cc +++ b/source/level_set_okz_reinitialization.cc @@ -258,10 +258,11 @@ struct ReinitializationMatrix template void -LevelSetOKZSolverReinitialization::reinitialize(const double dt, - const unsigned int stab_steps, - const unsigned int diff_steps, - const bool) +LevelSetOKZSolverReinitialization::reinitialize( + const double dt, + const unsigned int stab_steps, + const unsigned int diff_steps, + const std::function &compute_normal) { this->time_stepping.set_time_step(dt); @@ -292,7 +293,7 @@ LevelSetOKZSolverReinitialization::reinitialize(const double dt, { first_reinit_step = (tau == actual_diff_steps); if (first_reinit_step) - normal_operator.compute_normal(true); + compute_normal(true); // compute right hand side VectorType &rhs = this->system_rhs; From 20e23b294cfb1ea367f09f24e7c7ac863cd82291 Mon Sep 17 00:00:00 2001 From: Magdalena Schreter <70260016+mschreter@users.noreply.github.com> Date: Wed, 16 Dec 2020 12:27:00 +0100 Subject: [PATCH 15/16] Make reinitialize not virtual --- include/adaflo/level_set_okz_reinitialization.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/adaflo/level_set_okz_reinitialization.h b/include/adaflo/level_set_okz_reinitialization.h index a862e839..b4aed17b 100644 --- a/include/adaflo/level_set_okz_reinitialization.h +++ b/include/adaflo/level_set_okz_reinitialization.h @@ -100,7 +100,7 @@ class LevelSetOKZSolverReinitialization {} // performs reinitialization - virtual void + void reinitialize(const double dt, const unsigned int stab_steps, const unsigned int diff_steps = 0, From fe4f8e0c36640741c32f69b4f926a3de31e0d52b Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 26 Dec 2020 20:48:10 +0100 Subject: [PATCH 16/16] Fix scoped timer --- source/level_set_okz.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/source/level_set_okz.cc b/source/level_set_okz.cc index 3ed6cba1..998fc3f5 100644 --- a/source/level_set_okz.cc +++ b/source/level_set_okz.cc @@ -464,7 +464,6 @@ LevelSetOKZSolver::compute_curvature(const bool diffuse_large_values) // compute the curvature directly from C, but that is less accurate. TODO: // include that variant by a parameter { - TimerOutput::Scope timer(*this->timer, "LS compute normal."); this->compute_normal(false); } {