diff --git a/CMakeLists.txt b/CMakeLists.txt index 957e0f90..8f52e407 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,7 +55,12 @@ DEAL_II_INITIALIZE_CACHED_VARIABLES() # Set the source files to be compiled SET( TARGET_SRC + source/level_set_okz_preconditioner.cc + source/level_set_okz_compute_curvature.cc + source/level_set_okz_advance_concentration.cc + source/level_set_okz_compute_normal.cc source/level_set_okz.cc + source/level_set_okz_reinitialization.cc source/phase_field.cc source/phase_field_local.cc source/navier_stokes_matrix.cc diff --git a/include/adaflo/flow_base_algorithm.h b/include/adaflo/flow_base_algorithm.h index 46c8481c..243b8d52 100644 --- a/include/adaflo/flow_base_algorithm.h +++ b/include/adaflo/flow_base_algorithm.h @@ -56,8 +56,7 @@ namespace helpers std::set symmetry; std::set no_slip; - std::set fluid_type_plus; - std::set fluid_type_minus; + std::map>> fluid_type; std::pair periodic_boundaries[dim]; }; diff --git a/include/adaflo/level_set_okz.h b/include/adaflo/level_set_okz.h index e4b0f95f..7a7d6ceb 100644 --- a/include/adaflo/level_set_okz.h +++ b/include/adaflo/level_set_okz.h @@ -19,7 +19,10 @@ #include #include - +#include +#include +#include +#include using namespace dealii; @@ -42,29 +45,6 @@ class LevelSetOKZSolver : public LevelSetBaseAlgorithm virtual ~LevelSetOKZSolver() {} - void - advance_concentration_vmult( - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src) const; - - void - compute_normal_vmult(LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &sr) const; - - void - compute_normal_vmult(LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &sr) const; - - void - compute_curvature_vmult(LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &srcc, - const bool apply_diffusion) const; - - void - reinitialization_vmult(LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const bool diffuse_only) const; - virtual void initialize_data_structures(); @@ -108,63 +88,6 @@ class LevelSetOKZSolver : public LevelSetBaseAlgorithm const LinearAlgebra::distributed::Vector &src, const std::pair & cell_range); - template - void - local_advance_concentration( - const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range) const; - - template - void - local_advance_concentration_rhs( - const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range); - template - void - local_compute_normal(const MatrixFree & data, - LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &src, - const std::pair &cell_range) const; - template - void - local_compute_normal_rhs(const MatrixFree & data, - LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const; - - // diffusion_setting: 0: both terms, 1: only mass, 2: only diffusion - template - void - local_compute_curvature(const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair &cell_range) const; - template - void - local_compute_curvature_rhs( - const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range) const; - - template - void - local_reinitialize(const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range) const; - - template - void - local_reinitialize_rhs(const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range); - void local_projection_matrix( const MatrixFree & data, @@ -179,20 +102,17 @@ class LevelSetOKZSolver : public LevelSetBaseAlgorithm std::shared_ptr> &scratch, const std::pair & cell_range); - AlignedVector> artificial_viscosities; - AlignedVector>> evaluated_convection; - bool first_reinit_step; - double global_max_velocity; - DiagonalPreconditioner preconditioner; - - // In case we can better combine float/double solvers at some point... - MatrixFree matrix_free_float; - AlignedVector> cell_diameters_float; - // GrowingVectorMemory > - // vectors_normal; DiagonalPreconditioner preconditioner_float; + bool first_reinit_step; + double global_max_velocity; + DiagonalPreconditioner preconditioner; std::shared_ptr projection_matrix; std::shared_ptr ilu_projection_matrix; + + std::unique_ptr> normal_operator; + std::unique_ptr> curvatur_operator; + std::unique_ptr> reinit_operator; + std::unique_ptr> advection_operator; }; diff --git a/include/adaflo/level_set_okz_advance_concentration.h b/include/adaflo/level_set_okz_advance_concentration.h new file mode 100644 index 00000000..886da539 --- /dev/null +++ b/include/adaflo/level_set_okz_advance_concentration.h @@ -0,0 +1,183 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#ifndef __adaflo_level_set_advance_concentration_h +#define __adaflo_level_set_advance_concentration_h + +#include + +#include + +#include +#include +#include +#include + +using namespace dealii; + +/** + * Parameters of the advection-concentration operator. + */ +struct LevelSetOKZSolverAdvanceConcentrationParameter +{ + /** + * TODO + */ + unsigned int dof_index_ls; + + /** + * TODO + */ + unsigned int dof_index_vel; + + /** + * TODO + */ + unsigned int quad_index; + + /** + * TODO + */ + bool convection_stabilization; + + /** + * TODO + */ + bool do_iteration; + + /** + * TODO + */ + double tol_nl_iteration; + + /** + * TODO + */ + TimeSteppingParameters time; +}; + +/** + * Boundary descriptors of the advection-concentration operator. + */ +template +struct LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor +{ + /** + * TODO + */ + std::map>> dirichlet; + + /** + * TODO + */ + std::set symmetry; +}; + +template +class LevelSetOKZSolverAdvanceConcentration +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + + LevelSetOKZSolverAdvanceConcentration( + VectorType & solution, + const VectorType & solution_old, + const VectorType & solution_old_old, + VectorType & increment, + VectorType & rhs, + 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, + const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor &boundary, + const MatrixFree & matrix_free, + const LevelSetOKZSolverAdvanceConcentrationParameter & parameters, + double & global_max_velocity, + const DiagonalPreconditioner &preconditioner); + + virtual void + advance_concentration(const double dt); + + void + advance_concentration_vmult(VectorType &dst, const VectorType &src) const; + +private: + template + void + local_advance_concentration( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const; + + template + void + local_advance_concentration_rhs( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range); + + /** + * Parameters + */ + const LevelSetOKZSolverAdvanceConcentrationParameter parameters; // [i] + + /** + * Vector section + */ + VectorType & solution; // [o] new ls solution + const VectorType &solution_old; // [i] old ls solution + const VectorType &solution_old_old; // [i] old ls solution + VectorType & increment; // [-] temp + VectorType & rhs; // [-] temp + + const VectorType &vel_solution; // [i] new velocity solution + const VectorType &vel_solution_old; // [i] old velocity solution + const VectorType &vel_solution_old_old; // [i] old velocity solution + + /** + * MatrixFree + */ + const MatrixFree & matrix_free; // [i] + const AffineConstraints &constraints; // [i] + + /** + * Utility + */ + const ConditionalOStream &pcout; // [i] + TimeStepping time_stepping; // [i] + + /** + * Physics section + */ + const double & global_omega_diameter; // [i] + const AlignedVector> &cell_diameters; // [i] + const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor boundary; // [i] + AlignedVector> artificial_viscosities; // [-] + double & global_max_velocity; // [o] + AlignedVector>> evaluated_convection; // [o] + + /** + * Solver section + */ + const DiagonalPreconditioner &preconditioner; // [i] preconditioner +}; + +#endif diff --git a/include/adaflo/level_set_okz_compute_curvature.h b/include/adaflo/level_set_okz_compute_curvature.h new file mode 100644 index 00000000..9916b36f --- /dev/null +++ b/include/adaflo/level_set_okz_compute_curvature.h @@ -0,0 +1,157 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#ifndef __adaflo_level_compute_curvature_h +#define __adaflo_level_compute_curvature_h + +#include + +#include + +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +/** + * Parameters of the advection-concentration operator. + */ +struct LevelSetOKZSolverComputeCurvatureParameter +{ + /** + * TODO + */ + unsigned int dof_index_ls; + + /** + * TODO + */ + unsigned int dof_index_curvature; + + /** + * TODO + */ + unsigned int dof_index_normal; + + /** + * TODO + */ + unsigned int quad_index; + + /** + * TODO + */ + double epsilon; + + /** + * TODO + */ + bool approximate_projections; + + /** + * TODO + */ + bool curvature_correction; +}; + +template +class LevelSetOKZSolverComputeCurvature +{ +public: + LevelSetOKZSolverComputeCurvature( + LevelSetOKZSolverComputeNormal & normal_operator, + const AlignedVector> & cell_diameters, + const LinearAlgebra::distributed::BlockVector &normal_vector_field, + const AffineConstraints & constraints_curvature, + const AffineConstraints & constraints, + const double & epsilon_used, + LinearAlgebra::distributed::Vector & system_rhs, + const LevelSetOKZSolverComputeCurvatureParameter & parameters, + LinearAlgebra::distributed::Vector & solution_curvature, + const LinearAlgebra::distributed::Vector & solution_ls, + const MatrixFree & matrix_free, + const DiagonalPreconditioner & preconditioner, + std::shared_ptr & projection_matrix, + std::shared_ptr & ilu_projection_matrix); + + virtual void + compute_curvature(const bool diffuse_large_values = false); + + void + compute_curvature_vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &srcc, + const bool apply_diffusion) const; + +private: + // diffusion_setting: 0: both terms, 1: only mass, 2: only diffusion + template + void + local_compute_curvature(const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + template + void + local_compute_curvature_rhs( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const; + + /** + * Parameters + */ + const LevelSetOKZSolverComputeCurvatureParameter parameters; // [i] + + /** + * Other operators. + */ + LevelSetOKZSolverComputeNormal &normal_operator; // [i] + + /** + * Vector section + */ + LinearAlgebra::distributed::Vector & solution_curvature; // [i] + LinearAlgebra::distributed::Vector & rhs; // [-] + const LinearAlgebra::distributed::Vector & solution_ls; // [i] + const LinearAlgebra::distributed::BlockVector &normal_vector_field; // [i] + + /** + * MatrixFree + */ + const MatrixFree & matrix_free; // [i] + const AffineConstraints &constraints_curvature; // [i] + const AffineConstraints &constraints; // [i] + + /** + * Physics section + */ + const AlignedVector> &cell_diameters; // [i] + const double & epsilon_used; // [i] + + /** + * Solver section + */ + const DiagonalPreconditioner & preconditioner; // [i] + std::shared_ptr &projection_matrix; // [i] + std::shared_ptr & ilu_projection_matrix; // [i] +}; + +#endif diff --git a/include/adaflo/level_set_okz_compute_normal.h b/include/adaflo/level_set_okz_compute_normal.h new file mode 100644 index 00000000..003a39ad --- /dev/null +++ b/include/adaflo/level_set_okz_compute_normal.h @@ -0,0 +1,138 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#ifndef __adaflo_level_set_compute_normal_h +#define __adaflo_level_set_compute_normal_h + +#include + +#include + +#include +#include +#include + + +using namespace dealii; + +/** + * Parameters of the advection-concentration operator. + */ +struct LevelSetOKZSolverComputeNormalParameter +{ + /** + * TODO + */ + unsigned int dof_index_ls; + + /** + * TODO + */ + unsigned int dof_index_normal; + + /** + * TODO + */ + unsigned int quad_index; + + /** + * TODO + */ + double epsilon; + + /** + * TODO + */ + + bool approximate_projections; +}; + +template +class LevelSetOKZSolverComputeNormal +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + using BlockVectorType = LinearAlgebra::distributed::BlockVector; + + LevelSetOKZSolverComputeNormal( + BlockVectorType & normal_vector_field, + BlockVectorType & normal_vector_rhs, + VectorType & solution, + const AlignedVector> & cell_diameters, + const double & epsilon_used, + const double & minimal_edge_length, + const AffineConstraints & constraints_normals, + const LevelSetOKZSolverComputeNormalParameter ¶meters, + const MatrixFree & matrix_free, + const DiagonalPreconditioner & preconditioner, + const std::shared_ptr & projection_matrix, + const std::shared_ptr & ilu_projection_matrix); + + virtual void + compute_normal(const bool fast_computation); + + void + compute_normal_vmult(BlockVectorType &dst, const BlockVectorType &sr) const; + +private: + template + void + local_compute_normal(const MatrixFree & data, + LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src, + const std::pair &cell_range) const; + + template + void + local_compute_normal_rhs(const MatrixFree & data, + BlockVectorType & dst, + const VectorType & src, + const std::pair &cell_range) const; + + /** + * Parameters + */ + const LevelSetOKZSolverComputeNormalParameter parameters; // [i] + + /** + * Vector section + */ + BlockVectorType & normal_vector_field; // [o] + BlockVectorType & normal_vector_rhs; // [-] + const VectorType &vel_solution; // [i] + + /** + * MatrixFree + */ + const MatrixFree & matrix_free; // [i] + const AffineConstraints &constraints_normals; // [i] + + /** + * Physics section + */ + const AlignedVector> &cell_diameters; // [i] + const double & epsilon_used; // [i] + const double & minimal_edge_length; // [i] + + /** + * Solver section + */ + const DiagonalPreconditioner & preconditioner; // [i] + const std::shared_ptr &projection_matrix; // [i] + const std::shared_ptr & ilu_projection_matrix; // [i] +}; + +#endif diff --git a/include/adaflo/level_set_okz_preconditioner.h b/include/adaflo/level_set_okz_preconditioner.h new file mode 100644 index 00000000..be05b866 --- /dev/null +++ b/include/adaflo/level_set_okz_preconditioner.h @@ -0,0 +1,117 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#ifndef __adaflo_level_set_okz_preconditioner_h +#define __adaflo_level_set_okz_preconditioner_h + +#include + +#include +#include + +using namespace dealii; + +template > +inline void +initialize_mass_matrix_diagonal( + const MatrixFree &matrix_free, + const AffineConstraints & hanging_node_constraints, + const unsigned int dof_index, + const unsigned int quad_index, + DiagonalPreconditioner & preconditioner) +{ + LinearAlgebra::distributed::Vector diagonal; + matrix_free.initialize_dof_vector(diagonal, dof_index); + + const auto &dof_handler = matrix_free.get_dof_handler(dof_index); + const auto &fe = dof_handler.get_fe(); + const auto &quadrature = matrix_free.get_quadrature(quad_index); + const auto &mapping = *matrix_free.get_mapping_info().mapping; + + { + diagonal = 0; + FEValues fe_values(mapping, fe, quadrature, update_values | update_JxW_values); + Vector local_rhs(fe.dofs_per_cell); + std::vector local_dof_indices(fe.dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + { + Number value = 0; + for (const auto q : fe_values.quadrature_point_indices()) + value += fe_values.shape_value(i, q) * fe_values.shape_value(i, q) * + fe_values.JxW(q); + local_rhs(i) = value; + } + cell->get_dof_indices(local_dof_indices); + hanging_node_constraints.distribute_local_to_global(local_rhs, + local_dof_indices, + diagonal); + } + diagonal.compress(VectorOperation::add); + preconditioner.reinit(diagonal); + } +} + + + +namespace AssemblyData +{ + struct Data + { + Data() + { + AssertThrow(false, ExcNotImplemented()); + } + + Data(const unsigned int size) + : matrices(VectorizedArray::size(), FullMatrix(size, size)) + , dof_indices(size) + {} + + Data(const Data &other) + : matrices(other.matrices) + , dof_indices(other.dof_indices) + {} + + std::vector> matrices; + std::vector dof_indices; + }; +} // namespace AssemblyData + + + +template > +void +initialize_projection_matrix( + const MatrixFree &matrix_free, + const AffineConstraints & constraints_normals, + const unsigned int dof_index, + const unsigned int quad_index, + const Number & epsilon_used, + const Number & epsilon, + const AlignedVector & cell_diameters, + BlockMatrixExtension & projection_matrix, + BlockILUExtension & ilu_projection_matrix); + +#endif diff --git a/include/adaflo/level_set_okz_reinitialization.h b/include/adaflo/level_set_okz_reinitialization.h new file mode 100644 index 00000000..686af095 --- /dev/null +++ b/include/adaflo/level_set_okz_reinitialization.h @@ -0,0 +1,175 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#ifndef __adaflo_level_set_reinitialization_h +#define __adaflo_level_set_reinitialization_h + +#include + +#include + +#include +#include +#include +#include +#include + +using namespace dealii; + +/** + * Parameters of the reinitialization operator. + */ +struct LevelSetOKZSolverReinitializationParameter +{ + /** + * TODO + */ + unsigned int dof_index_ls; + + /** + * TODO + */ + unsigned int dof_index_normal; + + /** + * TODO + */ + unsigned int quad_index; + + /** + * TODO + */ + bool do_iteration; + + /** + * TODO + */ + TimeSteppingParameters time; +}; + +template +class LevelSetOKZSolverReinitialization +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + using BlockVectorType = LinearAlgebra::distributed::BlockVector; + + LevelSetOKZSolverReinitialization( + LevelSetOKZSolverComputeNormal & normal_operator, + const BlockVectorType & normal_vector_field, + const AlignedVector> & cell_diameters, + const double & epsilon_used, + const double & minimal_edge_length, + const AffineConstraints & constraints, + VectorType & solution_update, + VectorType & solution, + VectorType & system_rhs, + const ConditionalOStream & pcout, + const DiagonalPreconditioner & preconditioner, + const std::pair & last_concentration_range, + const LevelSetOKZSolverReinitializationParameter ¶meters, + bool & first_reinit_step, + const MatrixFree & matrix_free) + : parameters(parameters) + , normal_operator(normal_operator) + , solution(solution) + , solution_update(solution_update) + , system_rhs(system_rhs) + , normal_vector_field(normal_vector_field) + , matrix_free(matrix_free) + , constraints(constraints) + , cell_diameters(cell_diameters) + , epsilon_used(epsilon_used) + , minimal_edge_length(minimal_edge_length) + , last_concentration_range(last_concentration_range) + , first_reinit_step(first_reinit_step) + , pcout(pcout) + , time_stepping(parameters.time) + , preconditioner(preconditioner) + {} + + // 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 + reinitialization_vmult(VectorType & dst, + const VectorType &src, + const bool diffuse_only) const; + +private: + template + void + local_reinitialize(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const; + + template + void + local_reinitialize_rhs(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range); + + /** + * Parameters + */ + const LevelSetOKZSolverReinitializationParameter parameters; + + /** + * Other operators. + */ + LevelSetOKZSolverComputeNormal &normal_operator; + + /** + * Vector section + */ + VectorType &solution; // [o] + VectorType &solution_update; // [-] + VectorType &system_rhs; // [-] + + const BlockVectorType &normal_vector_field; // [i]; + + /** + * MatrixFree + */ + const MatrixFree & matrix_free; // [i] + const AffineConstraints &constraints; // [i] + + const AlignedVector> & cell_diameters; // [i] + const double & epsilon_used; // [i] + const double & minimal_edge_length; // [i] + const std::pair & last_concentration_range; // [i] + bool & first_reinit_step; // [?] + AlignedVector>> evaluated_normal; // [-] + + /** + * Utility + */ + const ConditionalOStream &pcout; // [i] + TimeStepping time_stepping; // [-] + + /** + * Solver section + */ + const DiagonalPreconditioner &preconditioner; // [i] +}; + +#endif diff --git a/include/adaflo/level_set_okz_template_instantations.h b/include/adaflo/level_set_okz_template_instantations.h new file mode 100644 index 00000000..d80e9a58 --- /dev/null +++ b/include/adaflo/level_set_okz_template_instantations.h @@ -0,0 +1,71 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#ifndef __adaflo_level_set_okz_template_instantiations_h +#define __adaflo_level_set_okz_template_instantiations_h + +#define EXPAND_OPERATIONS(OPERATION) \ + const unsigned int degree_u = this->navier_stokes.get_dof_handler_u().get_fe().degree; \ + const unsigned int ls_degree = this->parameters.concentration_subdivisions; \ + \ + AssertThrow(degree_u >= 2 && degree_u <= 5, ExcNotImplemented()); \ + AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ + if (ls_degree == 1) \ + { \ + if (degree_u == 2) \ + OPERATION(1, 2); \ + else if (degree_u == 3) \ + OPERATION(1, 3); \ + else if (degree_u == 4) \ + OPERATION(1, 4); \ + else if (degree_u == 5) \ + OPERATION(1, 5); \ + } \ + else if (ls_degree == 2) \ + { \ + if (degree_u == 2) \ + OPERATION(2, 2); \ + else if (degree_u == 3) \ + OPERATION(2, 3); \ + else if (degree_u == 4) \ + OPERATION(2, 4); \ + else if (degree_u == 5) \ + OPERATION(2, 5); \ + } \ + else if (ls_degree == 3) \ + { \ + if (degree_u == 2) \ + OPERATION(3, 2); \ + else if (degree_u == 3) \ + OPERATION(3, 3); \ + else if (degree_u == 4) \ + OPERATION(3, 4); \ + else if (degree_u == 5) \ + OPERATION(3, 5); \ + } \ + else if (ls_degree == 4) \ + { \ + if (degree_u == 2) \ + OPERATION(4, 2); \ + else if (degree_u == 3) \ + OPERATION(4, 3); \ + else if (degree_u == 4) \ + OPERATION(4, 4); \ + else if (degree_u == 5) \ + OPERATION(4, 5); \ + } + +#endif diff --git a/include/adaflo/parameters.h b/include/adaflo/parameters.h index 31cdb4e1..24a786f0 100644 --- a/include/adaflo/parameters.h +++ b/include/adaflo/parameters.h @@ -96,16 +96,15 @@ struct FlowParameters unsigned int print_solution_fields; bool output_wall_times; - TimeStepping::Scheme time_step_scheme; - double start_time; - double end_time; - double time_step_size_start; - double time_stepping_cfl; - double time_stepping_coef2; - double time_step_tolerance; - double time_step_size_max; - double time_step_size_min; - + TimeSteppingParameters::Scheme time_step_scheme; + double start_time; + double end_time; + double time_step_size_start; + double time_stepping_cfl; + double time_stepping_coef2; + double time_step_tolerance; + double time_step_size_max; + double time_step_size_min; // Two-phase specific parameters double density_diff; diff --git a/include/adaflo/time_stepping.h b/include/adaflo/time_stepping.h index b57bb404..5a8ce5af 100644 --- a/include/adaflo/time_stepping.h +++ b/include/adaflo/time_stepping.h @@ -28,26 +28,42 @@ using namespace dealii; struct FlowParameters; -class TimeStepping : public Subscriptor +struct TimeSteppingParameters { -public: - enum Scheme + /** + * TODO + */ + enum class Scheme { implicit_euler, explicit_euler, crank_nicolson, bdf_2 }; + Scheme time_step_scheme; + double start_time; + double end_time; + double time_step_size_start; + double time_stepping_cfl; + double time_stepping_coef2; + double time_step_tolerance; + double time_step_size_max; + double time_step_size_min; +}; + +class TimeStepping : public Subscriptor +{ +public: TimeStepping(const FlowParameters ¶meters); + TimeStepping(const TimeSteppingParameters ¶meters); + double start() const; double final() const; double - tolerance() const; - double step_size() const; double max_step_size() const; @@ -98,11 +114,9 @@ class TimeStepping : public Subscriptor set_final_time(double); void set_time_step(double); - void - set_tolerance(double); std::string name() const; - Scheme + TimeSteppingParameters::Scheme scheme() const; void @@ -111,32 +125,63 @@ class TimeStepping : public Subscriptor at_end() const; private: - double start_val; - double final_val; - double tolerance_val; - Scheme scheme_val; - double start_step_val; - double max_step_val; - double min_step_val; - double current_step_val; - double last_step_val; - double step_val; - double weight_val; - double weight_old_val; - double weight_old_old_val; - double factor_extrapol_old; - double factor_extrapol_old_old; - unsigned int step_no_val; - bool at_end_val; - bool weight_changed; - - double now_val; - double prev_val; - double tau1_val; - double tau2_val; + double start_val; // [m] start time; may be modified by set_start_time @todo: + // modification is nowhere used, keep? + double final_val; // [m] end time; may be modified by set_final_time @todo: modification + // is nowhere used, keep? + const TimeSteppingParameters::Scheme scheme_val; // [i] time integration scheme + const double start_step_val; // [i] initial value of the time increment + const double max_step_val; // [i] maximum value of the time increment + const double min_step_val; // [i] minimum value of the time increment + double current_step_val; // [m] current value of the time increment; + // initialized in the constructor by start_step_val + // can be modified in set_time_step(desired_value) + // fulfilling the criteria + // - 0.5 * step_size_prev <= current_step_val <= + // 2*step_size_prev + // - min_step_val <= current_step_val <= max_step_val + double last_step_val; // [m] constructor and restart() sets this parameter to zero. + // next() sets this parameter equal to current_step_val + // (corresponding to the previous time increment) + double step_val; // [m] current value of the time increment; m] + // - initialized in the constructor by start_step_val + // - changed in set_time_step to be equal to current_step_val + // @todo - ambiguous with current_step_val (?) + double weight_val; // [m] 1/time_increment + // - constructor: 1/start_step_val; + // - next(): 1/current_step_val + double weight_old_val; // [m] old time increment; + // - BDF2: this parameter is used for the time integration + // - else this parameter is used to determine weight_changed + double weight_old_old_val; // [m] old old time increment; only used in case of BDF2 + double factor_extrapol_old; // [m] extrapolation factor determined between the current + // and the last + // value of the increment + double factor_extrapol_old_old; // [m] extrapolation factor determined by the ratio + // between the current and + // and the old value of the time increment + unsigned int step_no_val; // [m] - constructor/restart(): initialize to zero + // - incremented by 1 in next() + bool at_end_val; // [m] determines if the end time is reached + bool weight_changed; // [m] determines if the integration weight has changed; this + // parameter is never reused + // this is used in navier_stokes.cc or phase_field.cc + double now_val; // [m] current time; time to be reached after time integration (t_n) + // - constructor/restart(): initialize to start_val + // - calculated as return argument from next() + double prev_val; // [m] old time; time at the begin of the integration step + double tau1_val; // [i] integration weight for multiplication with the current function + // value, i.e. f(t_n) + // - constructor: parameter depends on the scheme_val + double tau2_val; // [i] integration weight for multiplication with the old function + // value, i.e. f(t_n-1) + // - constructor: parameter depends on the scheme_val }; +/** + * Getter functions + */ inline double TimeStepping::start() const { @@ -151,7 +196,6 @@ TimeStepping::final() const } - inline double TimeStepping::step_size() const { @@ -173,13 +217,6 @@ TimeStepping::old_step_size() const return last_step_val; } -inline double -TimeStepping::tolerance() const -{ - return tolerance_val; -} - - inline double TimeStepping::now() const { @@ -219,7 +256,7 @@ TimeStepping::weight() const inline double TimeStepping::max_weight_uniform() const { - if (scheme_val == bdf_2) + if (scheme_val == TimeSteppingParameters::Scheme::bdf_2) return 1.5 / current_step_val; else return 1. / current_step_val; @@ -263,14 +300,6 @@ TimeStepping::set_final_time(double t) final_val = t; } - -inline void -TimeStepping::set_tolerance(double t) -{ - tolerance_val = t; -} - - inline bool TimeStepping::at_end() const { diff --git a/source/flow_base_algorithm.cc b/source/flow_base_algorithm.cc index 85f299d8..52a68c96 100644 --- a/source/flow_base_algorithm.cc +++ b/source/flow_base_algorithm.cc @@ -86,10 +86,12 @@ FlowBaseAlgorithm::set_velocity_dirichlet_boundary( case 0: break; case 1: - boundary->fluid_type_plus.insert(boundary_id); + boundary->fluid_type[boundary_id] = + std::make_shared>(1, 1); break; case -1: - boundary->fluid_type_minus.insert(boundary_id); + boundary->fluid_type[boundary_id] = + std::make_shared>(-1, 1); break; default: AssertThrow(false, ExcMessage("Unknown fluid type")); @@ -120,10 +122,12 @@ FlowBaseAlgorithm::set_open_boundary( case 0: break; case 1: - boundary->fluid_type_plus.insert(boundary_id); + boundary->fluid_type[boundary_id] = + std::make_shared>(1, 1); break; case -1: - boundary->fluid_type_minus.insert(boundary_id); + boundary->fluid_type[boundary_id] = + std::make_shared>(-1, 1); break; default: AssertThrow(false, ExcMessage("Unknown fluid type")); @@ -155,10 +159,12 @@ FlowBaseAlgorithm::set_open_boundary_with_normal_flux( case 0: break; case 1: - boundary->fluid_type_plus.insert(boundary_id); + boundary->fluid_type[boundary_id] = + std::make_shared>(1, 1); break; case -1: - boundary->fluid_type_minus.insert(boundary_id); + boundary->fluid_type[boundary_id] = + std::make_shared>(-1, 1); break; default: AssertThrow(false, ExcMessage("Unknown fluid type")); diff --git a/source/level_set_base.cc b/source/level_set_base.cc index ff4f6243..e533809b 100644 --- a/source/level_set_base.cc +++ b/source/level_set_base.cc @@ -121,16 +121,8 @@ LevelSetBaseAlgorithm::initialize_data_structures() // conditions on open boundaries Functions::ZeroFunction zero_func(1); std::map *> homogeneous_dirichlet; - for (typename std::set::const_iterator it = - this->boundary->fluid_type_plus.begin(); - it != this->boundary->fluid_type_plus.end(); - ++it) - homogeneous_dirichlet[*it] = &zero_func; - for (typename std::set::const_iterator it = - this->boundary->fluid_type_minus.begin(); - it != this->boundary->fluid_type_minus.end(); - ++it) - homogeneous_dirichlet[*it] = &zero_func; + for (const auto &it : this->boundary->fluid_type) + homogeneous_dirichlet[it.first] = &zero_func; VectorTools::interpolate_boundary_values(this->dof_handler, homogeneous_dirichlet, this->constraints); diff --git a/source/level_set_okz.cc b/source/level_set_okz.cc index 53d96ad8..e3ba337c 100644 --- a/source/level_set_okz.cc +++ b/source/level_set_okz.cc @@ -42,6 +42,8 @@ #include #include +#include +#include #include #include @@ -57,7 +59,146 @@ LevelSetOKZSolver::LevelSetOKZSolver(const FlowParameters ¶meters_in, Triangulation & tria_in) : LevelSetBaseAlgorithm(parameters_in, tria_in) , first_reinit_step(true) -{} +{ + { + LevelSetOKZSolverComputeNormalParameter params; + params.dof_index_ls = 2; + params.dof_index_normal = 4; + params.quad_index = 2; + params.epsilon = this->parameters.epsilon; + params.approximate_projections = this->parameters.approximate_projections; + + this->normal_operator = + std::make_unique>(this->normal_vector_field, + this->normal_vector_rhs, + this->solution.block(0), + this->cell_diameters, + this->epsilon_used, + this->minimal_edge_length, + this->constraints_normals, + params, + this->matrix_free, + this->preconditioner, + this->projection_matrix, + this->ilu_projection_matrix); + } + + { + LevelSetOKZSolverReinitializationParameter params; + + params.dof_index_ls = 2; + params.dof_index_normal = 4; + params.quad_index = 2; + params.do_iteration = this->parameters.do_iteration; + + // set time stepping parameters of level set to correspond with the values from + // Navier-Stokes + // @todo + params.time.time_step_scheme = this->parameters.time_step_scheme; + params.time.start_time = this->parameters.start_time; + params.time.end_time = this->parameters.end_time; + params.time.time_step_size_start = this->parameters.time_step_size_start; + params.time.time_stepping_cfl = this->parameters.time_stepping_cfl; + params.time.time_stepping_coef2 = this->parameters.time_stepping_coef2; + params.time.time_step_tolerance = this->parameters.time_step_tolerance; + params.time.time_step_size_max = this->parameters.time_step_size_max; + 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, + this->minimal_edge_length, + this->constraints, + this->solution_update.block(0), + this->solution.block(0), + this->system_rhs.block(0), + this->pcout, + this->preconditioner, + this->last_concentration_range, + params, + this->first_reinit_step, + this->matrix_free); + } + + { + LevelSetOKZSolverComputeCurvatureParameter params; + params.dof_index_ls = 2; + params.dof_index_curvature = 3; + params.dof_index_normal = 4; + params.quad_index = 2; + params.epsilon = this->parameters.epsilon; + params.approximate_projections = this->parameters.approximate_projections; + 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, + this->constraints, + this->epsilon_used, + this->system_rhs.block(0), + params, + this->solution.block(1), + this->solution.block(0), + this->matrix_free, + preconditioner, + projection_matrix, + ilu_projection_matrix); + } + + // set up advection operator + { + LevelSetOKZSolverAdvanceConcentrationParameter params; + + params.dof_index_ls = 2; + params.dof_index_vel = 0; + params.quad_index = 2; + params.convection_stabilization = this->parameters.convection_stabilization; + params.do_iteration = this->parameters.do_iteration; + params.tol_nl_iteration = this->parameters.tol_nl_iteration; + + LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor bcs; + + bcs.dirichlet = this->boundary->fluid_type; + bcs.symmetry = this->boundary->symmetry; + + // set time stepping parameters of level set to correspond with the values from + // Navier-Stokes + // @todo + params.time.time_step_scheme = this->parameters.time_step_scheme; + params.time.start_time = this->parameters.start_time; + params.time.end_time = this->parameters.end_time; + params.time.time_step_size_start = this->parameters.time_step_size_start; + params.time.time_stepping_cfl = this->parameters.time_stepping_cfl; + params.time.time_stepping_coef2 = this->parameters.time_stepping_coef2; + params.time.time_step_tolerance = this->parameters.time_step_tolerance; + params.time.time_step_size_max = this->parameters.time_step_size_max; + params.time.time_step_size_min = this->parameters.time_step_size_min; + + this->advection_operator = + std::make_unique>( + this->solution.block(0), + this->solution_old.block(0), + this->solution_old_old.block(0), + this->solution_update.block(0), + this->system_rhs.block(0), + 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); + } +} @@ -73,172 +214,28 @@ LevelSetOKZSolver::transform_distance_function( } - -namespace AssemblyData -{ - struct Data - { - Data() - { - AssertThrow(false, ExcNotImplemented()); - } - - Data(const unsigned int size) - : matrices(VectorizedArray::size(), FullMatrix(size, size)) - , dof_indices(size) - {} - - Data(const Data &other) - : matrices(other.matrices) - , dof_indices(other.dof_indices) - {} - - std::vector> matrices; - std::vector dof_indices; - }; -} // namespace AssemblyData - - // @sect4{LevelSetOKZSolver::make_grid_and_dofs} template void LevelSetOKZSolver::initialize_data_structures() { - cell_diameters_float.clear(); - matrix_free_float.clear(); this->LevelSetBaseAlgorithm::initialize_data_structures(); - artificial_viscosities.resize(this->matrix_free.n_cell_batches()); - evaluated_convection.resize(this->matrix_free.n_cell_batches() * - this->matrix_free.get_n_q_points(2)); - - // create diagonal preconditioner vector by assembly of mass matrix diagonal - LinearAlgebra::distributed::Vector diagonal(this->solution_update.block(0)); - { - diagonal = 0; - QIterated quadrature(QGauss<1>(2), this->parameters.concentration_subdivisions); - FEValues fe_values(this->mapping, - *this->fe, - quadrature, - update_values | update_JxW_values); - Vector local_rhs(this->fe->dofs_per_cell); - std::vector local_dof_indices(this->fe->dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = - this->dof_handler.begin_active(), - end = this->dof_handler.end(); - for (; cell != end; ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit(cell); - for (unsigned int i = 0; i < this->fe->dofs_per_cell; ++i) - { - double value = 0; - for (unsigned int q = 0; q < quadrature.size(); ++q) - value += fe_values.shape_value(i, q) * fe_values.shape_value(i, q) * - fe_values.JxW(q); - local_rhs(i) = value; - } - cell->get_dof_indices(local_dof_indices); - this->hanging_node_constraints.distribute_local_to_global(local_rhs, - local_dof_indices, - diagonal); - } - diagonal.compress(VectorOperation::add); - preconditioner.reinit(diagonal); - } - - // Initialize float matrix-free object for normal computation in case we - // want to do that at some point... - - // vectors_normal.release_unused_memory(); - // typename MatrixFree::AdditionalData data; - // data.tasks_parallel_scheme = - // MatrixFree::AdditionalData::partition_partition; - // data.mapping_update_flags = update_JxW_values | update_gradients; - // data.store_plain_indices = false; - // data.mpi_communicator = this->triangulation.get_communicator(); - // matrix_free_float.reinit(this->mapping, this->dof_handler, - // this->constraints, - // QIterated<1> (QGauss<1>(2), - // this->parameters.concentration_subdivisions), - // data); - // cell_diameters_float.resize(matrix_free_float.n_cell_batches()); - // std::map,double> diameters; - // for (unsigned int c=0; cmatrix_free.n_cell_batches(); ++c) - // for (unsigned int v=0; vmatrix_free.n_active_entries_per_cell_batch(c); ++v) - // diameters[std::make_pair(this->matrix_free.get_cell_iterator(c,v)->level(), - // this->matrix_free.get_cell_iterator(c,v)->index())] - // = this->cell_diameters[c][v]; - // for (unsigned int c=0; clevel(), - // matrix_free_float.get_cell_iterator(c,v)->index())]; - - // LinearAlgebra::distributed::Vector diagonal_f; - // diagonal_f = this->solution_update.block(0); - // preconditioner_float.reinit(diagonal_f); - - - // create sparse matrix for projection systems. - // - // First off is the creation of a mask that only adds those entries of - // FE_Q_iso_Q0 that are going to have a non-zero matrix entry -> this - // ensures as compact a matrix as for Q1 on the fine mesh. To find them, - // check terms in a mass matrix. - Table<2, bool> dof_mask(this->fe->dofs_per_cell, this->fe->dofs_per_cell); - { - QIterated quadrature(QGauss<1>(1), this->parameters.concentration_subdivisions); - FEValues fe_values(this->mapping, *this->fe, quadrature, update_values); - fe_values.reinit(this->dof_handler.begin()); - for (unsigned int i = 0; i < this->fe->dofs_per_cell; ++i) - for (unsigned int j = 0; j < this->fe->dofs_per_cell; ++j) - { - double sum = 0; - for (unsigned int q = 0; q < quadrature.size(); ++q) - sum += fe_values.shape_value(i, q) * fe_values.shape_value(j, q); - if (sum != 0) - dof_mask(i, j) = true; - } - } - { - IndexSet relevant_dofs; - DoFTools::extract_locally_relevant_dofs(this->dof_handler, relevant_dofs); - TrilinosWrappers::SparsityPattern csp; - csp.reinit(this->dof_handler.locally_owned_dofs(), - this->dof_handler.locally_owned_dofs(), - relevant_dofs, - get_communicator(this->triangulation)); - std::vector local_dof_indices(this->fe->dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = - this->dof_handler.begin_active(), - endc = this->dof_handler.end(); - for (; cell != endc; ++cell) - if (cell->is_locally_owned()) - { - cell->get_dof_indices(local_dof_indices); - this->constraints_normals.add_entries_local_to_global(local_dof_indices, - csp, - false, - dof_mask); - } - csp.compress(); - projection_matrix = std::make_shared(); - projection_matrix->reinit(csp); - } - { - AssemblyData::Data scratch_data(this->fe->dofs_per_cell); - auto scratch_local = - std::make_shared>(scratch_data); - unsigned int dummy = 0; - this->matrix_free.cell_loop(&LevelSetOKZSolver::local_projection_matrix, - this, - scratch_local, - dummy); - projection_matrix->compress(VectorOperation::add); - ilu_projection_matrix = std::make_shared(); - ilu_projection_matrix->initialize(*projection_matrix); - } + initialize_mass_matrix_diagonal( + this->matrix_free, this->hanging_node_constraints, 2, 2, preconditioner); + + projection_matrix = std::make_shared(); + ilu_projection_matrix = std::make_shared(); + + initialize_projection_matrix(this->matrix_free, + this->constraints_normals, + 2, + 2, + this->epsilon_used, + this->parameters.epsilon, + this->cell_diameters, + *projection_matrix, + *ilu_projection_matrix); } @@ -412,482 +409,6 @@ LevelSetOKZSolver::local_compute_force( } } - - -template -template -void -LevelSetOKZSolver::local_advance_concentration( - const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range) const -{ - // The second input argument below refers to which constrains should be used, - // 2 means constraints (for LS-function) - FEEvaluation ls_values(data, 2, 2); - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - const Tensor<1, dim, VectorizedArray> *velocities = - &evaluated_convection[cell * ls_values.n_q_points]; - ls_values.reinit(cell); - - ls_values.read_dof_values(src); - ls_values.evaluate(true, true); - - for (unsigned int q = 0; q < ls_values.n_q_points; ++q) - { - const VectorizedArray ls_val = ls_values.get_value(q); - const Tensor<1, dim, VectorizedArray> ls_grad = - ls_values.get_gradient(q); - ls_values.submit_value(ls_val * this->time_stepping.weight() + - ls_grad * velocities[q], - q); - if (this->parameters.convection_stabilization) - ls_values.submit_gradient(artificial_viscosities[cell] * ls_grad, q); - } - ls_values.integrate(true, this->parameters.convection_stabilization); - ls_values.distribute_local_to_global(dst); - } -} - - - -template -template -void -LevelSetOKZSolver::local_advance_concentration_rhs( - const MatrixFree & data, - LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &, - const std::pair &cell_range) -{ - // The second input argument below refers to which constrains should be used, - // 2 means constraints (for LS-function) and 0 means - // &navier_stokes.get_constraints_u() - FEEvaluation ls_values(data, 2, 2); - FEEvaluation ls_values_old(data, 2, 2); - FEEvaluation ls_values_old_old(data, 2, 2); - FEEvaluation vel_values(data, 0, 2); - FEEvaluation vel_values_old(data, 0, 2); - FEEvaluation vel_values_old_old(data, 0, 2); - - typedef VectorizedArray vector_t; - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - ls_values.reinit(cell); - ls_values_old.reinit(cell); - ls_values_old_old.reinit(cell); - vel_values.reinit(cell); - vel_values_old.reinit(cell); - vel_values_old_old.reinit(cell); - - vel_values.read_dof_values_plain(this->navier_stokes.solution.block(0)); - vel_values_old.read_dof_values_plain(this->navier_stokes.solution_old.block(0)); - vel_values_old_old.read_dof_values_plain( - this->navier_stokes.solution_old_old.block(0)); - ls_values.read_dof_values_plain(this->solution.block(0)); - ls_values_old.read_dof_values_plain(this->solution_old.block(0)); - ls_values_old_old.read_dof_values_plain(this->solution_old_old.block(0)); - - vel_values.evaluate(true, false); - vel_values_old.evaluate(true, false); - vel_values_old_old.evaluate(true, false); - ls_values.evaluate(true, true); - ls_values_old.evaluate(true, true); - ls_values_old_old.evaluate(true, true); - - if (this->parameters.convection_stabilization) - { - vector_t max_residual = vector_t(), max_velocity = vector_t(); - for (unsigned int q = 0; q < ls_values.n_q_points; ++q) - { - // compute residual of concentration equation - Tensor<1, dim, vector_t> u = - (vel_values_old.get_value(q) + vel_values_old_old.get_value(q)); - vector_t dc_dt = - (ls_values_old.get_value(q) - ls_values_old_old.get_value(q)) / - this->time_stepping.old_step_size(); - vector_t residual = std::abs( - dc_dt + - u * (ls_values_old.get_gradient(q) + ls_values_old_old.get_gradient(q)) * - 0.25); - max_residual = std::max(residual, max_residual); - max_velocity = std::max(std::sqrt(u * u), max_velocity); - } - double global_scaling = global_max_velocity * 2 * this->global_omega_diameter; - const vector_t cell_diameter = this->cell_diameters[cell]; - - artificial_viscosities[cell] = - 0.03 * max_velocity * cell_diameter * - std::min(make_vectorized_array(1.), 1. * max_residual / global_scaling); - } - - Tensor<1, dim, VectorizedArray> *velocities = - &evaluated_convection[cell * ls_values.n_q_points]; - - for (unsigned int q = 0; q < ls_values.n_q_points; ++q) - { - // compute right hand side - vector_t old_value = - this->time_stepping.weight_old() * ls_values_old.get_value(q); - if (this->time_stepping.scheme() == TimeStepping::bdf_2 && - this->time_stepping.step_no() > 1) - old_value += - this->time_stepping.weight_old_old() * ls_values_old_old.get_value(q); - const vector_t ls_val = ls_values.get_value(q); - const Tensor<1, dim, vector_t> ls_grad = ls_values.get_gradient(q); - vector_t residual = -(ls_val * this->time_stepping.weight() + - vel_values.get_value(q) * ls_grad + old_value); - ls_values.submit_value(residual, q); - if (this->parameters.convection_stabilization) - ls_values.submit_gradient(-artificial_viscosities[cell] * ls_grad, q); - velocities[q] = vel_values.get_value(q); - } - ls_values.integrate(true, this->parameters.convection_stabilization); - ls_values.distribute_local_to_global(dst); - } -} - - - -template -template -void -LevelSetOKZSolver::local_compute_normal( - const MatrixFree & data, - LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &src, - const std::pair & cell_range) const -{ - bool do_float = std::is_same::value; - // The second input argument below refers to which constrains should be used, - // 4 means constraints_normals - FEEvaluation phi(data, - do_float ? 0 : 4, - do_float ? 0 : 2); - 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 - const VectorizedArray *cell_diameters = - do_float ? - reinterpret_cast *>(cell_diameters_float.begin()) : - reinterpret_cast *>(this->cell_diameters.begin()); - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - phi.reinit(cell); - phi.read_dof_values(src); - phi.evaluate(true, true); - const VectorizedArray damping = - Number(4.) * - 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) - { - phi.submit_value(phi.get_value(q), q); - phi.submit_gradient(phi.get_gradient(q) * damping, q); - } - phi.integrate(true, true); - phi.distribute_local_to_global(dst); - } -} - - - -template -template -void -LevelSetOKZSolver::local_compute_normal_rhs( - const MatrixFree & data, - LinearAlgebra::distributed::BlockVector &dst, - const LinearAlgebra::distributed::Vector &, - const std::pair &cell_range) const -{ - // The second input argument below refers to which constrains should be used, - // 4 means constraints_normals and 2 means constraints (for LS-function) - FEEvaluation normal_values(data, 4, 2); - FEEvaluation ls_values(data, 2, 2); - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - normal_values.reinit(cell); - ls_values.reinit(cell); - - ls_values.read_dof_values_plain(this->solution.block(0)); - ls_values.evaluate(false, true, false); - - for (unsigned int q = 0; q < normal_values.n_q_points; ++q) - normal_values.submit_value(ls_values.get_gradient(q), q); - - normal_values.integrate(true, false); - normal_values.distribute_local_to_global(dst); - } -} - - - -template -template -void -LevelSetOKZSolver::local_compute_curvature( - const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range) const -{ - // The second input argument below refers to which constrains should be used, - // 3 means constraints_curvature - FEEvaluation phi(data, 3, 2); - const VectorizedArray min_diameter = - make_vectorized_array(this->epsilon_used / this->parameters.epsilon); - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - phi.reinit(cell); - phi.read_dof_values(src); - // If diffusion_setting is true a damping term is added to the weak form - // i.e. diffusion_setting=1 => diffusion_setting%2 == 1 is true. - phi.evaluate(diffusion_setting < 2, diffusion_setting % 2 == 1); - const VectorizedArray damping = - diffusion_setting % 2 == 1 ? - Utilities::fixed_power<2>( - std::max(min_diameter, - this->cell_diameters[cell] / static_cast(ls_degree))) : - VectorizedArray(); - for (unsigned int q = 0; q < phi.n_q_points; ++q) - { - if (diffusion_setting < 2) - phi.submit_value(phi.get_value(q), q); - if (diffusion_setting % 2 == 1) - phi.submit_gradient(phi.get_gradient(q) * damping, q); - } - phi.integrate(diffusion_setting < 2, diffusion_setting % 2 == 1); - phi.distribute_local_to_global(dst); - } -} - - - -template -template -void -LevelSetOKZSolver::local_compute_curvature_rhs( - const MatrixFree & data, - LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &, - const std::pair &cell_range) const -{ - // The second input argument below refers to which constrains should be used, - // 4 means constraints_normals and 3 constraints_curvature - FEEvaluation normal_values(data, 4, 2); - FEEvaluation curv_values(data, 3, 2); - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - normal_values.reinit(cell); - curv_values.reinit(cell); - - normal_values.read_dof_values_plain(this->normal_vector_field); - - // This computes (v, \nabla \cdot (n/|n|)). Since we store only \nabla - // \phi in the vector normal_values, the normalization must be done here - // and be done before we apply the derivatives, which is done in the - // code below. - 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); - } - - 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.integrate(true, false); - curv_values.distribute_local_to_global(dst); - } - } -} - - - -template -template -void -LevelSetOKZSolver::local_reinitialize( - const MatrixFree & data, - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const std::pair & cell_range) const -{ - const double dtau_inv = std::max(0.95 / (1. / (dim * dim) * this->minimal_edge_length / - this->parameters.concentration_subdivisions), - 1. / (5. * this->time_stepping.step_size())); - - // The second input argument below refers to which constrains should be used, - // 2 means constraints (for LS-function) - FEEvaluation phi(data, 2, 2); - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - phi.reinit(cell); - phi.read_dof_values(src); - phi.evaluate(true, true, false); - - VectorizedArray cell_diameter = this->cell_diameters[cell]; - VectorizedArray diffusion = - std::max(make_vectorized_array(this->epsilon_used), - cell_diameter / static_cast(ls_degree)); - - const Tensor<1, dim, VectorizedArray> *normal = - &evaluated_convection[cell * phi.n_q_points]; - for (unsigned int q = 0; q < phi.n_q_points; ++q) - if (!diffuse_only) - { - phi.submit_value(dtau_inv * phi.get_value(q), q); - phi.submit_gradient((diffusion * (normal[q] * phi.get_gradient(q))) * - normal[q], - q); - } - else - { - phi.submit_value(dtau_inv * phi.get_value(q), q); - phi.submit_gradient(phi.get_gradient(q) * diffusion, q); - } - - phi.integrate(true, true); - phi.distribute_local_to_global(dst); - } -} - - - -template -template -void -LevelSetOKZSolver::local_reinitialize_rhs( - const MatrixFree & data, - LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &, - const std::pair &cell_range) -{ - // The second input argument below refers to which constrains should be used, - // 2 means constraints (for LS-function) and 4 means constraints_normals - FEEvaluation phi(data, 2, 2); - FEEvaluation normals(data, 4, 2); - - for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) - { - phi.reinit(cell); - phi.read_dof_values_plain(this->solution.block(0)); - phi.evaluate(true, true, false); - - normals.reinit(cell); - normals.read_dof_values_plain(this->normal_vector_field); - normals.evaluate(true, false, false); - - VectorizedArray cell_diameter = this->cell_diameters[cell]; - VectorizedArray diffusion = - std::max(make_vectorized_array(this->epsilon_used), - cell_diameter / static_cast(ls_degree)); - - for (unsigned int q = 0; q < phi.n_q_points; ++q) - if (!diffuse_only) - { - 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()); - evaluated_convection[cell * phi.n_q_points + q] = normal; - } - // take normal as it was for the first reinit step - Tensor<1, dim, VectorizedArray> normal = - evaluated_convection[cell * phi.n_q_points + q]; - phi.submit_gradient(normal * - (0.5 * (1. - phi.get_value(q) * phi.get_value(q)) - - (normal * grad * diffusion)), - q); - } - else - { - phi.submit_gradient(-diffusion * phi.get_gradient(q), q); - } - - phi.integrate(false, true); - phi.distribute_local_to_global(dst); - } -} - - - -#define EXPAND_OPERATIONS(OPERATION) \ - const unsigned int degree_u = this->navier_stokes.get_dof_handler_u().get_fe().degree; \ - const unsigned int ls_degree = this->parameters.concentration_subdivisions; \ - \ - AssertThrow(degree_u >= 2 && degree_u <= 5, ExcNotImplemented()); \ - AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ - if (ls_degree == 1) \ - { \ - if (degree_u == 2) \ - OPERATION(1, 2); \ - else if (degree_u == 3) \ - OPERATION(1, 3); \ - else if (degree_u == 4) \ - OPERATION(1, 4); \ - else if (degree_u == 5) \ - OPERATION(1, 5); \ - } \ - else if (ls_degree == 2) \ - { \ - if (degree_u == 2) \ - OPERATION(2, 2); \ - else if (degree_u == 3) \ - OPERATION(2, 3); \ - else if (degree_u == 4) \ - OPERATION(2, 4); \ - else if (degree_u == 5) \ - OPERATION(2, 5); \ - } \ - else if (ls_degree == 3) \ - { \ - if (degree_u == 2) \ - OPERATION(3, 2); \ - else if (degree_u == 3) \ - OPERATION(3, 3); \ - else if (degree_u == 4) \ - OPERATION(3, 4); \ - else if (degree_u == 5) \ - OPERATION(3, 5); \ - } \ - else if (ls_degree == 4) \ - { \ - if (degree_u == 2) \ - OPERATION(4, 2); \ - else if (degree_u == 3) \ - OPERATION(4, 3); \ - else if (degree_u == 4) \ - OPERATION(4, 4); \ - else if (degree_u == 5) \ - OPERATION(4, 5); \ - } - - template void LevelSetOKZSolver::compute_force() @@ -913,510 +434,13 @@ LevelSetOKZSolver::compute_force() -template -void -LevelSetOKZSolver::advance_concentration_vmult( - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src) const -{ - dst = 0.; -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_advance_concentration, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - - - if (this->parameters.convection_stabilization) - { - // Boundary part of stabilization-term: - FEFaceValues fe_face_values(*this->fe, - this->matrix_free.get_face_quadrature(1), - update_values | update_gradients | - update_JxW_values | update_normal_vectors); - Vector cell_rhs(this->fe->dofs_per_cell); - std::vector local_dof_indices(this->fe->dofs_per_cell); - std::vector> local_gradients(fe_face_values.get_quadrature().size()); - src.update_ghost_values(); - - for (unsigned int mcell = 0; mcell < this->matrix_free.n_cell_batches(); ++mcell) - for (unsigned int v = 0; - v < this->matrix_free.n_active_entries_per_cell_batch(mcell); - ++v) - { - typename DoFHandler::active_cell_iterator cell = - this->matrix_free.get_cell_iterator(mcell, v, 2); - cell_rhs = 0; - - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) - { - if (cell->face(face)->at_boundary()) - { - if (this->boundary->symmetry.find(cell->face(face)->boundary_id()) == - this->boundary->symmetry.end()) - { - fe_face_values.reinit(cell, face); - fe_face_values.get_function_gradients(src, local_gradients); - for (unsigned int i = 0; i < this->fe->dofs_per_cell; ++i) - { - for (unsigned int q = 0; - q < fe_face_values.get_quadrature().size(); - ++q) - { - cell_rhs(i) += -((fe_face_values.shape_value(i, q) * - fe_face_values.normal_vector(q) * - artificial_viscosities[mcell][v] * - local_gradients[q]) * - fe_face_values.JxW(q)); - } - } - } - } - } - - cell->get_dof_indices(local_dof_indices); - this->constraints.distribute_local_to_global(cell_rhs, - local_dof_indices, - dst); - } - - dst.compress(VectorOperation::add); - } - - for (unsigned int i = 0; i < this->matrix_free.get_constrained_dofs(2).size(); ++i) - dst.local_element(this->matrix_free.get_constrained_dofs(2)[i]) = - preconditioner.get_vector().local_element( - this->matrix_free.get_constrained_dofs(2)[i]) * - src.local_element(this->matrix_free.get_constrained_dofs(2)[i]); -} - - - -template -struct AdvanceConcentrationMatrix -{ - AdvanceConcentrationMatrix(const LevelSetOKZSolver &problem) - : problem(problem) - {} - - void - vmult(LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src) const - { - problem.advance_concentration_vmult(dst, src); - } - - const LevelSetOKZSolver &problem; -}; - - - -template -void -LevelSetOKZSolver::compute_normal_vmult( - LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &src) const -{ - dst = 0.; -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_compute_normal, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - - // The number "4" below is so that constraints_normals is used - for (unsigned int i = 0; i < this->matrix_free.get_constrained_dofs(4).size(); ++i) - for (unsigned int d = 0; d < dim; ++d) - dst.block(d).local_element(this->matrix_free.get_constrained_dofs(4)[i]) = - preconditioner.get_vector().local_element( - this->matrix_free.get_constrained_dofs(4)[i]) * - src.block(d).local_element(this->matrix_free.get_constrained_dofs(4)[i]); -} - - - -template -void -LevelSetOKZSolver::compute_normal_vmult( - LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &src) const -{ - dst = 0.; -#define OPERATION(c_degree, u_degree) \ - matrix_free_float.cell_loop( \ - &LevelSetOKZSolver::template local_compute_normal, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - - // The number "4" below is so that constraints_normals is used - for (unsigned int i = 0; i < this->matrix_free.get_constrained_dofs(4).size(); ++i) - for (unsigned int d = 0; d < dim; ++d) - dst.block(d).local_element(this->matrix_free.get_constrained_dofs(4)[i]) = - preconditioner.get_vector().local_element( - this->matrix_free.get_constrained_dofs(4)[i]) * - src.block(d).local_element(this->matrix_free.get_constrained_dofs(4)[i]); -} - - - -template -struct ComputeNormalMatrix -{ - ComputeNormalMatrix(const LevelSetOKZSolver &problem) - : problem(problem) - {} - - template - void - vmult(LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &src) const - { - problem.compute_normal_vmult(dst, src); - } - - const LevelSetOKZSolver &problem; -}; - - - -template -class InverseNormalMatrix -{ -public: - InverseNormalMatrix( - GrowingVectorMemory> &mem, - const ComputeNormalMatrix & matrix, - const DiagonalPreconditioner & preconditioner) - : memory(mem) - , matrix(matrix) - , preconditioner(preconditioner) - {} - - void - vmult(LinearAlgebra::distributed::BlockVector & dst, - const LinearAlgebra::distributed::BlockVector &src) const - { - LinearAlgebra::distributed::BlockVector *src_f = memory.alloc(); - LinearAlgebra::distributed::BlockVector *dst_f = memory.alloc(); - - src_f->reinit(src); - dst_f->reinit(dst); - - *dst_f = 0; - *src_f = src; - ReductionControl control(10000, 1e-30, 1e-1); - SolverCG> solver(control, memory); - try - { - solver.solve(matrix, *dst_f, *src_f, preconditioner); - } - catch (...) - { - std::cout << "Error, normal solver did not converge!" << std::endl; - } - dst = *dst_f; - - memory.free(src_f); - memory.free(dst_f); - } - -private: - GrowingVectorMemory> &memory; - const ComputeNormalMatrix & matrix; - const DiagonalPreconditioner & preconditioner; -}; - - -template -void -LevelSetOKZSolver::compute_curvature_vmult( - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const bool apply_diffusion) const -{ - dst = 0.; - if (apply_diffusion) - { - // diffusion_setting will be 1 (true) in local_compute_curvature so that - // damping will be added -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_compute_curvature, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - } - else - { - // diffusion_setting will be 0 (fals) in local_compute_curvature so that - // NO damping will be added -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_compute_curvature, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - } - - // The numer "3" below is so that constraints_curvature is used - for (unsigned int i = 0; i < this->matrix_free.get_constrained_dofs(3).size(); ++i) - dst.local_element(this->matrix_free.get_constrained_dofs(3)[i]) = - preconditioner.get_vector().local_element( - this->matrix_free.get_constrained_dofs(3)[i]) * - src.local_element(this->matrix_free.get_constrained_dofs(3)[i]); -} - - - -template -struct ComputeCurvatureMatrix -{ - ComputeCurvatureMatrix(const LevelSetOKZSolver &problem) - : problem(problem) - {} - - void - vmult(LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src) const - { - problem.compute_curvature_vmult(dst, src, true); - } - - const LevelSetOKZSolver &problem; -}; - - - -template -void -LevelSetOKZSolver::reinitialization_vmult( - LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src, - const bool diffuse_only) const -{ - dst = 0.; - if (diffuse_only) - { -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_reinitialize, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - } - else - { -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_reinitialize, \ - this, \ - dst, \ - src) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - } - - for (unsigned int i = 0; i < this->matrix_free.get_constrained_dofs(2).size(); ++i) - dst.local_element(this->matrix_free.get_constrained_dofs(2)[i]) = - preconditioner.get_vector().local_element( - this->matrix_free.get_constrained_dofs(2)[i]) * - src.local_element(this->matrix_free.get_constrained_dofs(2)[i]); -} - - - -template -struct ReinitializationMatrix -{ - ReinitializationMatrix(const LevelSetOKZSolver &problem, const bool diffuse_only) - : problem(problem) - , diffuse_only(diffuse_only) - {} - - void - vmult(LinearAlgebra::distributed::Vector & dst, - const LinearAlgebra::distributed::Vector &src) const - { - problem.reinitialization_vmult(dst, src, diffuse_only); - } - - const LevelSetOKZSolver &problem; - const bool diffuse_only; -}; - - - // @sect4{LevelSetOKZSolver::advance_concentration} template void LevelSetOKZSolver::advance_concentration() { TimerOutput::Scope timer(*this->timer, "LS advance concentration."); - - // apply boundary values - { - std::map *> dirichlet; - Functions::ConstantFunction plus_func(1., 1); - for (typename std::set::const_iterator it = - this->boundary->fluid_type_plus.begin(); - it != this->boundary->fluid_type_plus.end(); - ++it) - dirichlet[*it] = &plus_func; - Functions::ConstantFunction minus_func(-1., 1); - for (typename std::set::const_iterator it = - this->boundary->fluid_type_minus.begin(); - it != this->boundary->fluid_type_minus.end(); - ++it) - dirichlet[*it] = &minus_func; - - std::map boundary_values; - VectorTools::interpolate_boundary_values(this->mapping, - this->dof_handler, - dirichlet, - boundary_values); - - for (typename std::map::const_iterator it = - boundary_values.begin(); - it != boundary_values.end(); - ++it) - if (this->solution.block(0).locally_owned_elements().is_element(it->first)) - this->solution.block(0)(it->first) = it->second; - this->solution.block(0).update_ghost_values(); - } - - // compute right hand side - global_max_velocity = this->get_maximal_velocity(); - LinearAlgebra::distributed::Vector &rhs = this->system_rhs.block(0); - LinearAlgebra::distributed::Vector &increment = this->solution_update.block(0); - rhs = 0; - -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_advance_concentration_rhs, \ - this, \ - rhs, \ - this->solution.block(0)) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - - AdvanceConcentrationMatrix matrix(*this); - - - - if (this->parameters.convection_stabilization) - { - // Boundary part of stabilization-term: - FEFaceValues fe_face_values(this->mapping, - *this->fe, - this->matrix_free.get_face_quadrature(1), - update_values | update_gradients | - update_JxW_values | update_normal_vectors); - - Vector cell_rhs(this->fe->dofs_per_cell); - std::vector local_dof_indices(this->fe->dofs_per_cell); - std::vector> local_gradients(fe_face_values.get_quadrature().size()); - - for (unsigned int mcell = 0; mcell < this->matrix_free.n_cell_batches(); ++mcell) - for (unsigned int v = 0; - v < this->matrix_free.n_active_entries_per_cell_batch(mcell); - ++v) - { - typename DoFHandler::active_cell_iterator cell = - this->matrix_free.get_cell_iterator(mcell, v, 2); - cell_rhs = 0; - - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) - { - if (cell->face(face)->at_boundary()) - { - if (this->boundary->symmetry.find(cell->face(face)->boundary_id()) == - this->boundary->symmetry.end()) - { - fe_face_values.reinit(cell, face); - fe_face_values.get_function_gradients(this->solution.block(0), - local_gradients); - - for (unsigned int i = 0; i < this->fe->dofs_per_cell; ++i) - { - for (unsigned int q = 0; - q < fe_face_values.get_quadrature().size(); - ++q) - { - cell_rhs(i) += ((fe_face_values.shape_value(i, q) * - fe_face_values.normal_vector(q) * - artificial_viscosities[mcell][v] * - local_gradients[q]) * - fe_face_values.JxW(q)); - } - } - } - } - } - - cell->get_dof_indices(local_dof_indices); - this->constraints.distribute_local_to_global(cell_rhs, - local_dof_indices, - this->system_rhs); - } - this->system_rhs.compress(VectorOperation::add); - } - - - // solve linear system with Bicgstab (non-symmetric system!) - unsigned int n_iterations = 0; - double initial_residual = 0.; - try - { - ReductionControl control(30, 0.05 * this->parameters.tol_nl_iteration, 1e-8); - SolverBicgstab>::AdditionalData - bicg_data; - bicg_data.exact_residual = false; - SolverBicgstab> solver(control, - bicg_data); - increment = 0; - solver.solve(matrix, increment, rhs, preconditioner); - n_iterations = control.last_step(); - initial_residual = control.initial_value(); - } - catch (const SolverControl::NoConvergence &) - { - // GMRES is typically slower but much more robust - ReductionControl control(3000, 0.05 * this->parameters.tol_nl_iteration, 1e-8); - SolverGMRES> solver(control); - solver.solve(matrix, increment, rhs, preconditioner); - n_iterations = 30 + control.last_step(); - } - if (!this->parameters.do_iteration) - this->pcout << " Concentration advance: advect [" << initial_residual << "/" - << n_iterations << "]"; - - this->constraints.distribute(increment); - this->solution.block(0) += increment; - this->solution.block(0).update_ghost_values(); + advection_operator->advance_concentration(this->time_stepping.step_size()); } @@ -1426,75 +450,8 @@ template void LevelSetOKZSolver::compute_normal(const bool fast_computation) { - // This function computes the normal from a projection of $\nabla C$ onto - // the space of linear finite elements (with some small damping) - TimerOutput::Scope timer(*this->timer, "LS compute normal."); - - // compute right hand side - this->normal_vector_rhs = 0; -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_compute_normal_rhs, \ - this, \ - this->normal_vector_rhs, \ - this->solution.block(0)) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - - if (this->parameters.approximate_projections == true) - { - // apply damping to avoid oscillations. this corresponds to one time - // step of exlicit Euler for a diffusion problem (need to avoid too - // large diffusions!) - const unsigned int n_steps = 3; - for (unsigned int i = 0; i < n_steps; ++i) - { - for (unsigned int block = 0; block < dim; ++block) - { - this->constraints_normals.distribute( - this->normal_vector_field.block(block)); - compute_curvature_vmult(this->normal_vector_rhs.block(block), - this->normal_vector_field.block(block), - 2); - } - preconditioner.vmult(this->normal_vector_rhs, this->normal_vector_rhs); - this->normal_vector_field.add(-0.05, this->normal_vector_rhs); - } - } - else - { - ComputeNormalMatrix matrix(*this); - - // solve linear system, reduce residual by 3e-7 in standard case and - // 1e-3 in fast case. Need a quite strict tolerance for the normal - // computation, otherwise the profile gets very bad when high curvatures - // appear and the solver fails - ReductionControl solver_control(4000, 1e-50, fast_computation ? 1e-5 : 1e-7); - - // ... in case we can somehow come up with a better combination of - // float/double solvers - // InverseNormalMatrix inverse(vectors_normal, matrix, - // preconditioner_float); - SolverCG> solver(solver_control); - // solver.solve (matrix, this->normal_vector_field, - // this->normal_vector_rhs, - // preconditioner); - solver.solve(*projection_matrix, - this->normal_vector_field, - this->normal_vector_rhs, - *ilu_projection_matrix); - // this->pcout << "N its normal: " << solver_control.last_step() << - // std::endl; - } - - - for (unsigned int d = 0; d < dim; ++d) - { - this->constraints_normals.distribute(this->normal_vector_field.block(d)); - this->normal_vector_field.block(d).update_ghost_values(); - } + normal_operator->compute_normal(fast_computation); } @@ -1502,76 +459,10 @@ LevelSetOKZSolver::compute_normal(const bool fast_computation) // @sect4{LevelSetOKZSolver::compute_normal} template void -LevelSetOKZSolver::compute_curvature(const bool) +LevelSetOKZSolver::compute_curvature(const bool 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 - compute_normal(false); - TimerOutput::Scope timer(*this->timer, "LS compute curvature."); - - // compute right hand side - LinearAlgebra::distributed::Vector &rhs = this->system_rhs.block(0); - rhs = 0; - -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_compute_curvature_rhs, \ - this, \ - rhs, \ - this->solution.block(0)) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - - // solve linear system for projection - if (this->parameters.approximate_projections == true) - preconditioner.vmult(this->solution.block(1), rhs); - else - { - ComputeCurvatureMatrix matrix(*this); - - ReductionControl solver_control(2000, 1e-50, 1e-8); - SolverCG> cg(solver_control); - // cg.solve (matrix, this->solution.block(1), rhs, preconditioner); - cg.solve(*projection_matrix, this->solution.block(1), rhs, *ilu_projection_matrix); - // this->pcout << "N its curv: " << solver_control.last_step() << - // std::endl; - } - - // correct curvature away from the zero level set by computing the distance - // and correcting the value, if so requested - if (this->parameters.curvature_correction == true) - { - for (unsigned int i = 0; i < this->solution.block(1).local_size(); ++i) - if (this->solution.block(1).local_element(i) > 1e-4) - { - const double c_val = this->solution.block(0).local_element(i); - const double distance = - (1 - c_val * c_val) > 1e-2 ? - this->epsilon_used * std::log((1. + c_val) / (1. - c_val)) : - 0; - - this->solution.block(1).local_element(i) = - 1. / (1. / this->solution.block(1).local_element(i) + distance / (dim - 1)); - } - } - - this->constraints_curvature.distribute(this->solution.block(1)); - - // apply damping to avoid oscillations. this corresponds to one time - // step of explicit Euler for a diffusion problem (need to avoid too - // large diffusions!) - if (this->parameters.approximate_projections == true) - for (unsigned int i = 0; i < 8; ++i) - { - compute_curvature_vmult(rhs, this->solution.block(1), 2); - preconditioner.vmult(rhs, rhs); - this->solution.block(1).add(-0.05, rhs); - this->constraints.distribute(this->solution.block(1)); - } - this->solution.block(1).update_ghost_values(); + curvatur_operator->compute_curvature(diffuse_large_values); } @@ -1647,100 +538,13 @@ template void LevelSetOKZSolver::reinitialize(const unsigned int stab_steps, const unsigned int diff_steps, - const bool) + const bool diffuse_cells_with_large_curvature_only) { - // This function assembles and solves for a given profile using the approach - // described in the paper by Olsson, Kreiss, and Zahedi. - - std::cout.precision(3); - - // perform several reinitialization steps until we reach the maximum number - // of steps. - // - // TODO: make an adaptive choice of the number of iterations - unsigned actual_diff_steps = diff_steps; - if (this->last_concentration_range.first < -1.02 || - this->last_concentration_range.second > 1.02) - actual_diff_steps += 3; - if (!this->parameters.do_iteration) - this->pcout << (this->time_stepping.now() == this->time_stepping.start() ? " " : - " and ") - << "reinitialize ("; - for (unsigned int tau = 0; tau < actual_diff_steps + stab_steps; tau++) - { - first_reinit_step = (tau == actual_diff_steps); - if (first_reinit_step) - compute_normal(true); - - TimerOutput::Scope timer(*this->timer, "LS reinitialization step."); - - // compute right hand side - LinearAlgebra::distributed::Vector &rhs = this->system_rhs.block(0); - LinearAlgebra::distributed::Vector &increment = - this->solution_update.block(0); - rhs = 0; - - if (tau < actual_diff_steps) - { -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_reinitialize_rhs, \ - this, \ - rhs, \ - this->solution.block(0)) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - } - else - { -#define OPERATION(c_degree, u_degree) \ - this->matrix_free.cell_loop( \ - &LevelSetOKZSolver::template local_reinitialize_rhs, \ - this, \ - rhs, \ - this->solution.block(0)) - - EXPAND_OPERATIONS(OPERATION); -#undef OPERATION - } - - // solve linear system - { - ReinitializationMatrix matrix(*this, tau < actual_diff_steps); - increment = 0; - - // reduce residual by 1e-6. To obtain good interface shapes, it is - // essential that this tolerance is relative to the rhs - // (ReductionControl steered solver, last argument determines the - // solver) - ReductionControl solver_control(2000, 1e-50, 1e-6); - SolverCG> cg(solver_control); - cg.solve(matrix, increment, rhs, preconditioner); - this->constraints.distribute(increment); - if (!this->parameters.do_iteration) - { - if (tau < actual_diff_steps) - this->pcout << "d" << solver_control.last_step(); - else - this->pcout << solver_control.last_step(); - } - } - - this->solution.block(0) += increment; - this->solution.block(0).update_ghost_values(); - - // check residual - const double update_norm = increment.l2_norm(); - if (update_norm < 1e-6) - break; - - if (!this->parameters.do_iteration && tau < actual_diff_steps + stab_steps - 1) - this->pcout << " + "; - } - - if (!this->parameters.do_iteration) - this->pcout << ")" << std::endl << std::flush; + 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); } diff --git a/source/level_set_okz_advance_concentration.cc b/source/level_set_okz_advance_concentration.cc new file mode 100644 index 00000000..4c862767 --- /dev/null +++ b/source/level_set_okz_advance_concentration.cc @@ -0,0 +1,593 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +namespace +{ + /** + * Compute maximal velocity for a given vector and the corresponding + * dof-handler object. + */ + template + double + get_maximal_velocity(const DoFHandler &dof_handler, + const VectorType & solution, + const Quadrature &quad_in) + { + // [PM] We use QIterated in the case of hex mesh for backwards compatibility. + const Quadrature quadrature_formula = + dof_handler.get_fe().reference_cell_type() == ReferenceCell::get_hypercube(dim) ? + Quadrature( + QIterated(QTrapez<1>(), dof_handler.get_fe().tensor_degree() + 1)) : + quad_in; + + FEValues fe_values(dof_handler.get_fe(), quadrature_formula, update_values); + std::vector> velocity_values(quadrature_formula.size()); + + const FEValuesExtractors::Vector velocities(0); + + double max_velocity = 0; + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + fe_values.reinit(cell); + fe_values[velocities].get_function_values(solution, velocity_values); + + for (const auto q : fe_values.quadrature_point_indices()) + max_velocity = std::max(max_velocity, velocity_values[q].norm()); + } + + return Utilities::MPI::max(max_velocity, get_communicator(dof_handler)); + } +} // namespace + + + +#define EXPAND_OPERATIONS(OPERATION) \ + if (this->matrix_free.get_dof_handler(parameters.dof_index_vel) \ + .get_fe() \ + .reference_cell_type() != ReferenceCell::get_hypercube(dim)) \ + { \ + OPERATION(-1, -1); \ + } \ + else \ + { \ + const unsigned int degree_u = \ + this->matrix_free.get_dof_handler(parameters.dof_index_vel) \ + .get_fe() \ + .tensor_degree(); \ + const unsigned int ls_degree = \ + this->matrix_free.get_dof_handler(parameters.dof_index_ls) \ + .get_fe() \ + .tensor_degree(); \ + \ + AssertThrow(degree_u >= 1 && degree_u <= 5, ExcNotImplemented()); \ + AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ + if (ls_degree == 1) \ + { \ + if (degree_u == 1) \ + OPERATION(1, 1); \ + else if (degree_u == 2) \ + OPERATION(1, 2); \ + else if (degree_u == 3) \ + OPERATION(1, 3); \ + else if (degree_u == 4) \ + OPERATION(1, 4); \ + else if (degree_u == 5) \ + OPERATION(1, 5); \ + } \ + else if (ls_degree == 2) \ + { \ + if (degree_u == 1) \ + OPERATION(2, 1); \ + else if (degree_u == 2) \ + OPERATION(2, 2); \ + else if (degree_u == 3) \ + OPERATION(2, 3); \ + else if (degree_u == 4) \ + OPERATION(2, 4); \ + else if (degree_u == 5) \ + OPERATION(2, 5); \ + } \ + else if (ls_degree == 3) \ + { \ + if (degree_u == 1) \ + OPERATION(3, 1); \ + else if (degree_u == 2) \ + OPERATION(3, 2); \ + else if (degree_u == 3) \ + OPERATION(3, 3); \ + else if (degree_u == 4) \ + OPERATION(3, 4); \ + else if (degree_u == 5) \ + OPERATION(3, 5); \ + } \ + else if (ls_degree == 4) \ + { \ + if (degree_u == 1) \ + OPERATION(4, 1); \ + else if (degree_u == 2) \ + OPERATION(4, 2); \ + else if (degree_u == 3) \ + OPERATION(4, 3); \ + else if (degree_u == 4) \ + OPERATION(4, 4); \ + else if (degree_u == 5) \ + OPERATION(4, 5); \ + } \ + } + + + +template +LevelSetOKZSolverAdvanceConcentration::LevelSetOKZSolverAdvanceConcentration( + VectorType & solution, + const VectorType & solution_old, + const VectorType & solution_old_old, + VectorType & increment, + VectorType & rhs, + 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, + const LevelSetOKZSolverAdvanceConcentrationBoundaryDescriptor &boundary, + const MatrixFree & matrix_free, + const LevelSetOKZSolverAdvanceConcentrationParameter & parameters, + double & global_max_velocity, + const DiagonalPreconditioner & preconditioner) + : parameters(parameters) + , solution(solution) + , solution_old(solution_old) + , solution_old_old(solution_old_old) + , increment(increment) + , rhs(rhs) + , vel_solution(vel_solution) + , vel_solution_old(vel_solution_old) + , vel_solution_old_old(vel_solution_old_old) + , matrix_free(matrix_free) + , constraints(constraints) + , pcout(pcout) + , time_stepping(parameters.time) + , global_omega_diameter(global_omega_diameter) + , cell_diameters(cell_diameters) + , boundary(boundary) + , global_max_velocity(global_max_velocity) + , preconditioner(preconditioner) +{} + + + +template +template +void +LevelSetOKZSolverAdvanceConcentration::local_advance_concentration( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const +{ + // The second input argument below refers to which constrains should be used, + // 2 means constraints (for LS-function) + const unsigned int n_q_points = ls_degree == -1 ? 0 : 2 * ls_degree; + + FEEvaluation ls_values(data, + parameters.dof_index_ls, + parameters.quad_index); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + const Tensor<1, dim, VectorizedArray> *velocities = + &evaluated_convection[cell * ls_values.n_q_points]; + ls_values.reinit(cell); + + ls_values.gather_evaluate(src, true, true); + + for (unsigned int q = 0; q < ls_values.n_q_points; ++q) + { + const VectorizedArray ls_val = ls_values.get_value(q); + const Tensor<1, dim, VectorizedArray> ls_grad = + ls_values.get_gradient(q); + ls_values.submit_value(ls_val * this->time_stepping.weight() + + ls_grad * velocities[q], + q); + if (this->parameters.convection_stabilization) + ls_values.submit_gradient(artificial_viscosities[cell] * ls_grad, q); + } + ls_values.integrate_scatter(true, this->parameters.convection_stabilization, dst); + } +} + + + +template +template +void +LevelSetOKZSolverAdvanceConcentration::local_advance_concentration_rhs( + const MatrixFree &data, + VectorType & dst, + const VectorType &, + const std::pair &cell_range) +{ + // The second input argument below refers to which constrains should be used, + // 2 means constraints (for LS-function) and 0 means + // &navier_stokes.get_constraints_u() + const unsigned int n_q_points = ls_degree == -1 ? 0 : 2 * ls_degree; + + FEEvaluation ls_values(data, + parameters.dof_index_ls, + parameters.quad_index); + FEEvaluation ls_values_old(data, + parameters.dof_index_ls, + parameters.quad_index); + FEEvaluation ls_values_old_old(data, + parameters.dof_index_ls, + parameters.quad_index); + FEEvaluation vel_values(data, + parameters.dof_index_vel, + parameters.quad_index); + FEEvaluation vel_values_old( + data, parameters.dof_index_vel, parameters.quad_index); + FEEvaluation vel_values_old_old( + data, parameters.dof_index_vel, parameters.quad_index); + + typedef VectorizedArray vector_t; + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + ls_values.reinit(cell); + ls_values_old.reinit(cell); + ls_values_old_old.reinit(cell); + vel_values.reinit(cell); + vel_values_old.reinit(cell); + vel_values_old_old.reinit(cell); + + vel_values.read_dof_values_plain(vel_solution); + vel_values_old.read_dof_values_plain(vel_solution_old); + vel_values_old_old.read_dof_values_plain(vel_solution_old_old); + ls_values.read_dof_values_plain(this->solution); + ls_values_old.read_dof_values_plain(this->solution_old); + ls_values_old_old.read_dof_values_plain(this->solution_old_old); + + vel_values.evaluate(true, false); + vel_values_old.evaluate(true, false); + vel_values_old_old.evaluate(true, false); + ls_values.evaluate(true, true); + ls_values_old.evaluate(true, true); + ls_values_old_old.evaluate(true, true); + + if (this->parameters.convection_stabilization) + { + vector_t max_residual = vector_t(), max_velocity = vector_t(); + for (unsigned int q = 0; q < ls_values.n_q_points; ++q) + { + // compute residual of concentration equation + Tensor<1, dim, vector_t> u = + (vel_values_old.get_value(q) + vel_values_old_old.get_value(q)); + vector_t dc_dt = + (ls_values_old.get_value(q) - ls_values_old_old.get_value(q)) / + this->time_stepping.old_step_size(); + vector_t residual = std::abs( + dc_dt + + u * (ls_values_old.get_gradient(q) + ls_values_old_old.get_gradient(q)) * + 0.25); + max_residual = std::max(residual, max_residual); + max_velocity = std::max(std::sqrt(u * u), max_velocity); + } + double global_scaling = global_max_velocity * 2 * this->global_omega_diameter; + const vector_t cell_diameter = this->cell_diameters[cell]; + + artificial_viscosities[cell] = + 0.03 * max_velocity * cell_diameter * + std::min(make_vectorized_array(1.), 1. * max_residual / global_scaling); + } + + Tensor<1, dim, VectorizedArray> *velocities = + &evaluated_convection[cell * ls_values.n_q_points]; + + for (unsigned int q = 0; q < ls_values.n_q_points; ++q) + { + // compute right hand side + vector_t old_value = + this->time_stepping.weight_old() * ls_values_old.get_value(q); + if (this->time_stepping.scheme() == TimeSteppingParameters::Scheme::bdf_2 && + this->time_stepping.step_no() > 1) + old_value += + this->time_stepping.weight_old_old() * ls_values_old_old.get_value(q); + const vector_t ls_val = ls_values.get_value(q); + const Tensor<1, dim, vector_t> ls_grad = ls_values.get_gradient(q); + vector_t residual = -(ls_val * this->time_stepping.weight() + + vel_values.get_value(q) * ls_grad + old_value); + ls_values.submit_value(residual, q); + if (this->parameters.convection_stabilization) + ls_values.submit_gradient(-artificial_viscosities[cell] * ls_grad, q); + velocities[q] = vel_values.get_value(q); + } + ls_values.integrate_scatter(true, this->parameters.convection_stabilization, dst); + } +} + + + +template +void +LevelSetOKZSolverAdvanceConcentration::advance_concentration_vmult( + VectorType & dst, + const VectorType &src) const +{ + dst = 0.; +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop( \ + &LevelSetOKZSolverAdvanceConcentration< \ + dim>::template local_advance_concentration, \ + this, \ + dst, \ + src) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + + + if (this->parameters.convection_stabilization) + { + const auto &dof_handler = + this->matrix_free.get_dof_handler(parameters.dof_index_ls); + const auto &fe = dof_handler.get_fe(); + + // Boundary part of stabilization-term: + FEFaceValues fe_face_values( + fe, + this->matrix_free.get_face_quadrature(parameters.quad_index), + update_values | update_gradients | update_JxW_values | update_normal_vectors); + Vector cell_rhs(fe.dofs_per_cell); + std::vector local_dof_indices(fe.dofs_per_cell); + std::vector> local_gradients(fe_face_values.get_quadrature().size()); + src.update_ghost_values(); + + for (unsigned int mcell = 0; mcell < this->matrix_free.n_cell_batches(); ++mcell) + for (unsigned int v = 0; + v < this->matrix_free.n_active_entries_per_cell_batch(mcell); + ++v) + { + typename DoFHandler::active_cell_iterator cell = + this->matrix_free.get_cell_iterator(mcell, + v, + this->parameters.dof_index_ls); + cell_rhs = 0; + + for (const auto &face : cell->face_iterators()) + { + if (face->at_boundary() == false) + continue; + + if (this->boundary.symmetry.find(face->boundary_id()) != + this->boundary.symmetry.end()) + continue; + + fe_face_values.reinit(cell, face); + fe_face_values.get_function_gradients(src, local_gradients); + for (const auto i : fe_face_values.dof_indices()) + for (const auto q : fe_face_values.quadrature_point_indices()) + cell_rhs(i) += + -((fe_face_values.shape_value(i, q) * + fe_face_values.normal_vector(q) * + artificial_viscosities[mcell][v] * local_gradients[q]) * + fe_face_values.JxW(q)); + } + + cell->get_dof_indices(local_dof_indices); + this->constraints.distribute_local_to_global(cell_rhs, + local_dof_indices, + dst); + } + + dst.compress(VectorOperation::add); + } + + for (unsigned int i = 0; + i < this->matrix_free.get_constrained_dofs(this->parameters.dof_index_ls).size(); + ++i) + dst.local_element( + this->matrix_free.get_constrained_dofs(this->parameters.dof_index_ls)[i]) = + preconditioner.get_vector().local_element( + this->matrix_free.get_constrained_dofs(this->parameters.dof_index_ls)[i]) * + src.local_element( + this->matrix_free.get_constrained_dofs(this->parameters.dof_index_ls)[i]); +} + + + +template +struct AdvanceConcentrationMatrix +{ + AdvanceConcentrationMatrix(const LevelSetOKZSolverAdvanceConcentration &problem) + : problem(problem) + {} + + void + vmult(VectorType &dst, const VectorType &src) const + { + problem.advance_concentration_vmult(dst, src); + } + + const LevelSetOKZSolverAdvanceConcentration &problem; +}; + + + +// @sect4{LevelSetOKZSolverAdvanceConcentration::advance_concentration} +template +void +LevelSetOKZSolverAdvanceConcentration::advance_concentration(const double dt) +{ + this->time_stepping.set_time_step(dt); + this->time_stepping.next(); + + if (evaluated_convection.size() != + this->matrix_free.n_cell_batches() * + this->matrix_free.get_n_q_points(parameters.quad_index)) + evaluated_convection.resize(this->matrix_free.n_cell_batches() * + this->matrix_free.get_n_q_points(parameters.quad_index)); + + const auto &mapping = *this->matrix_free.get_mapping_info().mapping; + const auto &dof_handler = this->matrix_free.get_dof_handler(parameters.dof_index_ls); + const auto &fe = dof_handler.get_fe(); + + if (artificial_viscosities.size() != this->matrix_free.n_cell_batches()) + artificial_viscosities.resize(this->matrix_free.n_cell_batches()); + + // apply boundary values + { + std::map *> dirichlet; + + for (const auto &b : this->boundary.dirichlet) + dirichlet[b.first] = b.second.get(); + + std::map boundary_values; + VectorTools::interpolate_boundary_values(mapping, + dof_handler, + dirichlet, + boundary_values); + + for (const auto &it : boundary_values) + if (this->solution.locally_owned_elements().is_element(it.first)) + this->solution(it.first) = it.second; + this->solution.update_ghost_values(); + } + + // compute right hand side + global_max_velocity = + get_maximal_velocity(matrix_free.get_dof_handler(parameters.dof_index_vel), + vel_solution, + matrix_free.get_quadrature(parameters.quad_index)); + rhs = 0; + +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop( \ + &LevelSetOKZSolverAdvanceConcentration< \ + dim>::template local_advance_concentration_rhs, \ + this, \ + rhs, \ + this->solution) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + + AdvanceConcentrationMatrix matrix(*this); + + + + if (this->parameters.convection_stabilization) + { + // Boundary part of stabilization-term: + FEFaceValues fe_face_values( + mapping, + fe, + this->matrix_free.get_face_quadrature(parameters.quad_index), + update_values | update_gradients | update_JxW_values | update_normal_vectors); + + Vector cell_rhs(fe.dofs_per_cell); + std::vector local_dof_indices(fe.dofs_per_cell); + std::vector> local_gradients(fe_face_values.get_quadrature().size()); + + for (unsigned int mcell = 0; mcell < this->matrix_free.n_cell_batches(); ++mcell) + for (unsigned int v = 0; + v < this->matrix_free.n_active_entries_per_cell_batch(mcell); + ++v) + { + typename DoFHandler::active_cell_iterator cell = + this->matrix_free.get_cell_iterator(mcell, + v, + this->parameters.dof_index_ls); + cell_rhs = 0; + + for (const auto face : cell->face_iterators()) + { + if (face->at_boundary() == false) + continue; + + if (this->boundary.symmetry.find(face->boundary_id()) != + this->boundary.symmetry.end()) + continue; + + fe_face_values.reinit(cell, face); + fe_face_values.get_function_gradients(this->solution, local_gradients); + + for (const auto i : fe_face_values.dof_indices()) + for (const auto q : fe_face_values.quadrature_point_indices()) + cell_rhs(i) += + ((fe_face_values.shape_value(i, q) * + fe_face_values.normal_vector(q) * + artificial_viscosities[mcell][v] * local_gradients[q]) * + fe_face_values.JxW(q)); + } + + cell->get_dof_indices(local_dof_indices); + this->constraints.distribute_local_to_global(cell_rhs, + local_dof_indices, + this->rhs); + } + this->rhs.compress(VectorOperation::add); + } + + + // solve linear system with Bicgstab (non-symmetric system!) + unsigned int n_iterations = 0; + double initial_residual = 0.; + try + { + ReductionControl control(30, 0.05 * this->parameters.tol_nl_iteration, 1e-8); + SolverBicgstab::AdditionalData bicg_data; + bicg_data.exact_residual = false; + SolverBicgstab solver(control, bicg_data); + increment = 0; + solver.solve(matrix, increment, rhs, preconditioner); + n_iterations = control.last_step(); + initial_residual = control.initial_value(); + } + catch (const SolverControl::NoConvergence &) + { + // GMRES is typically slower but much more robust + ReductionControl control(3000, 0.05 * this->parameters.tol_nl_iteration, 1e-8); + SolverGMRES solver(control); + solver.solve(matrix, increment, rhs, preconditioner); + n_iterations = 30 + control.last_step(); + } + if (!this->parameters.do_iteration) + this->pcout << " Concentration advance: advect [" << initial_residual << "/" + << n_iterations << "]"; + + this->constraints.distribute(increment); + this->solution += increment; + this->solution.update_ghost_values(); +} + + +template class LevelSetOKZSolverAdvanceConcentration<2>; +template class LevelSetOKZSolverAdvanceConcentration<3>; diff --git a/source/level_set_okz_compute_curvature.cc b/source/level_set_okz_compute_curvature.cc new file mode 100644 index 00000000..87f40fb6 --- /dev/null +++ b/source/level_set_okz_compute_curvature.cc @@ -0,0 +1,321 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + + +#include +#include + +#include + +#include + + +#define EXPAND_OPERATIONS(OPERATION) \ + const unsigned int ls_degree = \ + this->matrix_free.get_dof_handler(parameters.dof_index_ls).get_fe().tensor_degree(); \ + \ + AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ + if (ls_degree == 1) \ + OPERATION(1, 0); \ + else if (ls_degree == 2) \ + OPERATION(2, 0); \ + else if (ls_degree == 3) \ + OPERATION(3, 0); \ + else if (ls_degree == 4) \ + OPERATION(4, 0); + +template +LevelSetOKZSolverComputeCurvature::LevelSetOKZSolverComputeCurvature( + LevelSetOKZSolverComputeNormal & normal_operator, + const AlignedVector> & cell_diameters, + const LinearAlgebra::distributed::BlockVector &normal_vector_field, + const AffineConstraints & constraints_curvature, + const AffineConstraints & constraints, + const double & epsilon_used, + LinearAlgebra::distributed::Vector & system_rhs, + const LevelSetOKZSolverComputeCurvatureParameter & parameters, + LinearAlgebra::distributed::Vector & solution_curvature, + const LinearAlgebra::distributed::Vector & solution_ls, + const MatrixFree & matrix_free, + const DiagonalPreconditioner & preconditioner, + 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) + , normal_vector_field(normal_vector_field) + , matrix_free(matrix_free) + , constraints_curvature(constraints_curvature) + , constraints(constraints) + , cell_diameters(cell_diameters) + , epsilon_used(epsilon_used) + , preconditioner(preconditioner) + , projection_matrix(projection_matrix) + , ilu_projection_matrix(ilu_projection_matrix) +{} + + +template +template +void +LevelSetOKZSolverComputeCurvature::local_compute_curvature( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const +{ + // The second input argument below refers to which constrains should be used, + // 3 means constraints_curvature + FEEvaluation phi(data, + parameters.dof_index_curvature, + parameters.quad_index); + const VectorizedArray min_diameter = + make_vectorized_array(this->epsilon_used / this->parameters.epsilon); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + // If diffusion_setting is true a damping term is added to the weak form + // i.e. diffusion_setting=1 => diffusion_setting%2 == 1 is true. + phi.evaluate(diffusion_setting < 2, diffusion_setting % 2 == 1); + const VectorizedArray damping = + diffusion_setting % 2 == 1 ? + Utilities::fixed_power<2>( + std::max(min_diameter, + this->cell_diameters[cell] / static_cast(ls_degree))) : + VectorizedArray(); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + if (diffusion_setting < 2) + phi.submit_value(phi.get_value(q), q); + if (diffusion_setting % 2 == 1) + phi.submit_gradient(phi.get_gradient(q) * damping, q); + } + phi.integrate(diffusion_setting < 2, diffusion_setting % 2 == 1); + phi.distribute_local_to_global(dst); + } +} + + + +template +template +void +LevelSetOKZSolverComputeCurvature::local_compute_curvature_rhs( + const MatrixFree & data, + LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &, + const std::pair &cell_range) const +{ + // The second input argument below refers to which constrains should be used, + // 4 means constraints_normals and 3 constraints_curvature + FEEvaluation normal_values( + data, parameters.dof_index_normal, parameters.quad_index); + FEEvaluation curv_values( + data, parameters.dof_index_curvature, parameters.quad_index); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + normal_values.reinit(cell); + curv_values.reinit(cell); + + normal_values.read_dof_values_plain(this->normal_vector_field); + + // This computes (v, \nabla \cdot (n/|n|)). Since we store only \nabla + // \phi in the vector normal_values, the normalization must be done here + // and be done before we apply the derivatives, which is done in the + // code below. + 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); + } + + 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.integrate(true, false); + curv_values.distribute_local_to_global(dst); + } + } +} + + + +template +void +LevelSetOKZSolverComputeCurvature::compute_curvature_vmult( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const bool apply_diffusion) const +{ + dst = 0.; + if (apply_diffusion) + { + // diffusion_setting will be 1 (true) in local_compute_curvature so that + // damping will be added +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverComputeCurvature< \ + dim>::template local_compute_curvature, \ + this, \ + dst, \ + src) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + } + else + { + // diffusion_setting will be 0 (fals) in local_compute_curvature so that + // NO damping will be added +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverComputeCurvature< \ + dim>::template local_compute_curvature, \ + this, \ + dst, \ + src) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + } + + // The numer "3" below is so that constraints_curvature is used + for (unsigned int i = 0; + i < this->matrix_free.get_constrained_dofs(parameters.dof_index_curvature).size(); + ++i) + dst.local_element( + this->matrix_free.get_constrained_dofs(parameters.dof_index_curvature)[i]) = + preconditioner.get_vector().local_element( + this->matrix_free.get_constrained_dofs(parameters.dof_index_curvature)[i]) * + src.local_element( + this->matrix_free.get_constrained_dofs(parameters.dof_index_curvature)[i]); +} + + + +template +struct ComputeCurvatureMatrix +{ + ComputeCurvatureMatrix(const LevelSetOKZSolverComputeCurvature &problem) + : problem(problem) + {} + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + problem.compute_curvature_vmult(dst, src, true); + } + + const LevelSetOKZSolverComputeCurvature &problem; +}; + + + +// @sect4{LevelSetOKZSolverComputeCurvature::compute_normal} +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; + +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverComputeCurvature< \ + dim>::template local_compute_curvature_rhs, \ + this, \ + rhs, \ + this->solution_ls) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + + // solve linear system for projection + if (this->parameters.approximate_projections == true) + preconditioner.vmult(this->solution_curvature, rhs); + else + { + ComputeCurvatureMatrix matrix(*this); + + ReductionControl solver_control(2000, 1e-50, 1e-8); + SolverCG> cg(solver_control); + // cg.solve (matrix, this->solution_curvature, rhs, preconditioner); + cg.solve(*projection_matrix, this->solution_curvature, rhs, *ilu_projection_matrix); + // this->pcout << "N its curv: " << solver_control.last_step() << + // std::endl; + } + + // correct curvature away from the zero level set by computing the distance + // and correcting the value, if so requested + if (this->parameters.curvature_correction == true) + { + for (unsigned int i = 0; i < this->solution_curvature.local_size(); ++i) + if (this->solution_curvature.local_element(i) > 1e-4) + { + const double c_val = this->solution_ls.local_element(i); + const double distance = + (1 - c_val * c_val) > 1e-2 ? + this->epsilon_used * std::log((1. + c_val) / (1. - c_val)) : + 0; + + this->solution_curvature.local_element(i) = + 1. / + (1. / this->solution_curvature.local_element(i) + distance / (dim - 1)); + } + } + + this->constraints_curvature.distribute(this->solution_curvature); + + // apply damping to avoid oscillations. this corresponds to one time + // step of explicit Euler for a diffusion problem (need to avoid too + // large diffusions!) + if (this->parameters.approximate_projections == true) + for (unsigned int i = 0; i < 8; ++i) + { + compute_curvature_vmult(rhs, this->solution_curvature, 2); + preconditioner.vmult(rhs, rhs); + this->solution_curvature.add(-0.05, rhs); + this->constraints.distribute(this->solution_curvature); + } + this->solution_curvature.update_ghost_values(); +} + + + +template class LevelSetOKZSolverComputeCurvature<2>; +template class LevelSetOKZSolverComputeCurvature<3>; diff --git a/source/level_set_okz_compute_normal.cc b/source/level_set_okz_compute_normal.cc new file mode 100644 index 00000000..ff9f7872 --- /dev/null +++ b/source/level_set_okz_compute_normal.cc @@ -0,0 +1,261 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + + +#include +#include + +#include + +#include + + +#define EXPAND_OPERATIONS(OPERATION) \ + const unsigned int ls_degree = \ + this->matrix_free.get_dof_handler(parameters.dof_index_ls).get_fe().tensor_degree(); \ + \ + AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ + if (ls_degree == 1) \ + OPERATION(1, 0); \ + else if (ls_degree == 2) \ + OPERATION(2, 0); \ + else if (ls_degree == 3) \ + OPERATION(3, 0); \ + else if (ls_degree == 4) \ + OPERATION(4, 0); + + +template +LevelSetOKZSolverComputeNormal::LevelSetOKZSolverComputeNormal( + BlockVectorType & normal_vector_field, + BlockVectorType & normal_vector_rhs, + VectorType & solution, + const AlignedVector> & cell_diameters, + const double & epsilon_used, + const double & minimal_edge_length, + const AffineConstraints & constraints_normals, + const LevelSetOKZSolverComputeNormalParameter ¶meters, + const MatrixFree & matrix_free, + const DiagonalPreconditioner & preconditioner, + const std::shared_ptr & projection_matrix, + const std::shared_ptr & ilu_projection_matrix) + : parameters(parameters) + , normal_vector_field(normal_vector_field) + , normal_vector_rhs(normal_vector_rhs) + , vel_solution(solution) + , matrix_free(matrix_free) + , constraints_normals(constraints_normals) + , cell_diameters(cell_diameters) + , epsilon_used(epsilon_used) + , minimal_edge_length(minimal_edge_length) + , preconditioner(preconditioner) + , projection_matrix(projection_matrix) + , ilu_projection_matrix(ilu_projection_matrix) +{} + + +template +template +void +LevelSetOKZSolverComputeNormal::local_compute_normal( + const MatrixFree & data, + LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src, + const std::pair & cell_range) const +{ + // The second input argument below refers to which constrains should be used, + // 4 means constraints_normals + FEEvaluation phi(data, 4, 2); + 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 + const VectorizedArray *cell_diameters = this->cell_diameters.begin(); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(true, true); + const VectorizedArray damping = + Number(4.) * + 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) + { + phi.submit_value(phi.get_value(q), q); + phi.submit_gradient(phi.get_gradient(q) * damping, q); + } + phi.integrate(true, true); + phi.distribute_local_to_global(dst); + } +} + + + +template +template +void +LevelSetOKZSolverComputeNormal::local_compute_normal_rhs( + const MatrixFree &data, + BlockVectorType & dst, + const VectorType &, + const std::pair &cell_range) const +{ + // The second input argument below refers to which constrains should be used, + // 4 means constraints_normals and 2 means constraints (for LS-function) + FEEvaluation normal_values(data, 4, 2); + FEEvaluation ls_values(data, 2, 2); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + normal_values.reinit(cell); + ls_values.reinit(cell); + + ls_values.read_dof_values_plain(this->vel_solution); + ls_values.evaluate(false, true, false); + + for (unsigned int q = 0; q < normal_values.n_q_points; ++q) + normal_values.submit_value(ls_values.get_gradient(q), q); + + normal_values.integrate(true, false); + normal_values.distribute_local_to_global(dst); + } +} + + + +template +void +LevelSetOKZSolverComputeNormal::compute_normal_vmult( + BlockVectorType & dst, + const BlockVectorType &src) const +{ + dst = 0.; +#define OPERATION(c_degree, dummy) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverComputeNormal< \ + dim>::template local_compute_normal, \ + this, \ + dst, \ + src) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + + // The number "4" below is so that constraints_normals is used + for (unsigned int i = 0; i < this->matrix_free.get_constrained_dofs(4).size(); ++i) + for (unsigned int d = 0; d < dim; ++d) + dst.block(d).local_element(this->matrix_free.get_constrained_dofs(4)[i]) = + preconditioner.get_vector().local_element( + this->matrix_free.get_constrained_dofs(4)[i]) * + src.block(d).local_element(this->matrix_free.get_constrained_dofs(4)[i]); +} + + + +template +struct ComputeNormalMatrix +{ + ComputeNormalMatrix(const LevelSetOKZSolverComputeNormal &problem) + : problem(problem) + {} + + template + void + vmult(LinearAlgebra::distributed::BlockVector & dst, + const LinearAlgebra::distributed::BlockVector &src) const + { + problem.compute_normal_vmult(dst, src); + } + + const LevelSetOKZSolverComputeNormal &problem; +}; + + + +// @sect4{LevelSetOKZSolverComputeNormal::compute_normal} +template +void +LevelSetOKZSolverComputeNormal::compute_normal(const bool fast_computation) +{ + // This function computes the normal from a projection of $\nabla C$ onto + // the space of linear finite elements (with some small damping) + + // compute right hand side + this->normal_vector_rhs = 0; +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop( \ + &LevelSetOKZSolverComputeNormal::template local_compute_normal_rhs, \ + this, \ + this->normal_vector_rhs, \ + this->vel_solution) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + + if (this->parameters.approximate_projections == true) + { + AssertThrow(false, ExcNotImplemented()); + // [PM]: The following code has been removed since it was not used in any test + // and it introduced annoying cyclic dependencies. + + /* + // apply damping to avoid oscillations. this corresponds to one time + // step of exlicit Euler for a diffusion problem (need to avoid too + // large diffusions!) + const unsigned int n_steps = 3; + for (unsigned int i = 0; i < n_steps; ++i) + { + for (unsigned int block = 0; block < dim; ++block) + { + this->constraints_normals.distribute( + this->normal_vector_field.block(block)); + curvatur_operator.compute_curvature_vmult( + this->normal_vector_rhs.block(block), + this->normal_vector_field.block(block), + 2); + */ + } + else + { + ComputeNormalMatrix matrix(*this); + + // solve linear system, reduce residual by 3e-7 in standard case and + // 1e-3 in fast case. Need a quite strict tolerance for the normal + // computation, otherwise the profile gets very bad when high curvatures + // appear and the solver fails + ReductionControl solver_control(4000, 1e-50, fast_computation ? 1e-5 : 1e-7); + + SolverCG solver(solver_control); + // solver.solve (matrix, this->normal_vector_field, + // this->normal_vector_rhs, + // preconditioner); + solver.solve(*projection_matrix, + this->normal_vector_field, + this->normal_vector_rhs, + *ilu_projection_matrix); + // this->pcout << "N its normal: " << solver_control.last_step() << + // std::endl; + } + + + for (unsigned int d = 0; d < dim; ++d) + { + this->constraints_normals.distribute(this->normal_vector_field.block(d)); + this->normal_vector_field.block(d).update_ghost_values(); + } +} + +template class LevelSetOKZSolverComputeNormal<2>; +template class LevelSetOKZSolverComputeNormal<3>; diff --git a/source/level_set_okz_matrix.cc b/source/level_set_okz_matrix.cc index e8d742b3..cd94e6bf 100644 --- a/source/level_set_okz_matrix.cc +++ b/source/level_set_okz_matrix.cc @@ -377,18 +377,9 @@ LevelSetOKZMatrixSolver::advance_concentration() // apply boundary values { std::map *> dirichlet; - Functions::ConstantFunction plus_func(1., 1); - for (typename std::set::const_iterator it = - this->boundary->fluid_type_plus.begin(); - it != this->boundary->fluid_type_plus.end(); - ++it) - dirichlet[*it] = &plus_func; - Functions::ConstantFunction minus_func(-1., 1); - for (typename std::set::const_iterator it = - this->boundary->fluid_type_minus.begin(); - it != this->boundary->fluid_type_minus.end(); - ++it) - dirichlet[*it] = &minus_func; + + for (const auto &b : this->boundary->fluid_type) + dirichlet[b.first] = b.second.get(); std::map boundary_values; VectorTools::interpolate_boundary_values(this->dof_handler, @@ -502,7 +493,7 @@ LevelSetOKZMatrixSolver::advance_concentration() { const Tensor<1, dim> &velocity = velocity_values[q]; double old_value = -this->time_stepping.weight_old() * old_values[q]; - if (this->time_stepping.scheme() == TimeStepping::bdf_2 && + if (this->time_stepping.scheme() == TimeSteppingParameters::Scheme::bdf_2 && this->time_stepping.step_no() > 1) old_value -= this->time_stepping.weight_old_old() * old_old_values[q]; diff --git a/source/level_set_okz_preconditioner.cc b/source/level_set_okz_preconditioner.cc new file mode 100644 index 00000000..b9e986ce --- /dev/null +++ b/source/level_set_okz_preconditioner.cc @@ -0,0 +1,196 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include +#include + + + +#define EXPAND_OPERATIONS(OPERATION) \ + AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ + if (ls_degree == 1) \ + { \ + OPERATION(1, 0); \ + } \ + else if (ls_degree == 2) \ + { \ + OPERATION(2, 0); \ + } \ + else if (ls_degree == 3) \ + { \ + OPERATION(3, 0); \ + } \ + else if (ls_degree == 4) \ + { \ + OPERATION(4, 0); \ + } + +template +void +initialize_projection_matrix( + const MatrixFree &matrix_free, + const AffineConstraints & constraints_normals, + const unsigned int dof_index, + const unsigned int quad_index, + const Number & epsilon_used, + const Number & epsilon, + const AlignedVector & cell_diameters, + BlockMatrixExtension & projection_matrix, + BlockILUExtension & ilu_projection_matrix) +{ + const auto &dof_handler = matrix_free.get_dof_handler(dof_index); + const auto &fe = dof_handler.get_fe(); + const auto &quadrature = matrix_free.get_quadrature(quad_index); + const auto &mapping = *matrix_free.get_mapping_info().mapping; + + // create sparse matrix for projection systems. + // + // First off is the creation of a mask that only adds those entries of + // FE_Q_iso_Q0 that are going to have a non-zero matrix entry -> this + // ensures as compact a matrix as for Q1 on the fine mesh. To find them, + // check terms in a mass matrix. + Table<2, bool> dof_mask(fe.dofs_per_cell, fe.dofs_per_cell); + { + FEValues fe_values(mapping, fe, quadrature, update_values); + fe_values.reinit(dof_handler.begin()); + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) + { + double sum = 0; + for (unsigned int q = 0; q < quadrature.size(); ++q) + sum += fe_values.shape_value(i, q) * fe_values.shape_value(j, q); + if (sum != 0) + dof_mask(i, j) = true; + } + } + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, relevant_dofs); + TrilinosWrappers::SparsityPattern csp; + csp.reinit(dof_handler.locally_owned_dofs(), + dof_handler.locally_owned_dofs(), + relevant_dofs, + get_communicator(dof_handler)); + std::vector local_dof_indices(fe.dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + cell->get_dof_indices(local_dof_indices); + constraints_normals.add_entries_local_to_global(local_dof_indices, + csp, + false, + dof_mask); + } + csp.compress(); + projection_matrix.reinit(csp); + } + { + AssemblyData::Data scratch_data(fe.dofs_per_cell); + auto scratch_local = + std::make_shared>(scratch_data); + unsigned int dummy = 0; + matrix_free.template cell_loop< + std::shared_ptr>, + unsigned int>( + [&](const auto &data, auto &scratch_data, const auto &, const auto cell_range) { + const unsigned int ls_degree = fe.tensor_degree(); + +#define OPERATION(c_degree, u_degree) \ + FEEvaluation phi(data, dof_index, quad_index); \ + AssemblyData::Data & scratch = scratch_data->get(); \ + \ + const VectorizedArray min_diameter = \ + make_vectorized_array(epsilon_used / epsilon); \ + \ + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) \ + { \ + phi.reinit(cell); \ + const VectorizedArray damping = \ + 4. * \ + Utilities::fixed_power<2>( \ + std::max(min_diameter, \ + cell_diameters[cell] / static_cast(fe.tensor_degree()))); \ + \ + for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) \ + { \ + for (unsigned int j = 0; j < phi.dofs_per_cell; ++j) \ + phi.begin_dof_values()[j] = VectorizedArray(); \ + phi.begin_dof_values()[i] = 1.; \ + phi.evaluate(true, true); \ + for (unsigned int q = 0; q < phi.n_q_points; ++q) \ + { \ + phi.submit_value(phi.get_value(q), q); \ + phi.submit_gradient(phi.get_gradient(q) * damping, q); \ + } \ + phi.integrate(true, true); \ + for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell); ++v) \ + for (unsigned int j = 0; j < phi.dofs_per_cell; ++j) \ + scratch.matrices[v](phi.get_shape_info().lexicographic_numbering[j], \ + phi.get_shape_info().lexicographic_numbering[i]) = \ + phi.begin_dof_values()[j][v]; \ + } \ + for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell); ++v) \ + { \ + typename DoFHandler::active_cell_iterator dcell = \ + matrix_free.get_cell_iterator(cell, v, dof_index); \ + dcell->get_dof_indices(scratch.dof_indices); \ + constraints_normals.distribute_local_to_global( \ + scratch.matrices[v], \ + scratch.dof_indices, \ + static_cast(projection_matrix)); \ + } \ + } + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + }, + scratch_local, + dummy); + projection_matrix.compress(VectorOperation::add); + ilu_projection_matrix.initialize(projection_matrix); + } +} + +template void +initialize_projection_matrix<2, double, VectorizedArray>( + const MatrixFree<2, 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<3, double, VectorizedArray>( + const MatrixFree<3, 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); diff --git a/source/level_set_okz_reinitialization.cc b/source/level_set_okz_reinitialization.cc new file mode 100644 index 00000000..e329acd1 --- /dev/null +++ b/source/level_set_okz_reinitialization.cc @@ -0,0 +1,336 @@ +// -------------------------------------------------------------------------- +// +// Copyright (C) 2020 by the adaflo authors +// +// This file is part of the adaflo library. +// +// The adaflo library is free software; you can use it, redistribute it, +// and/or modify it under the terms of the GNU Lesser General Public License +// as published by the Free Software Foundation; either version 2.1 of the +// License, or (at your option) any later version. The full text of the +// license can be found in the file LICENSE at the top level of the adaflo +// distribution. +// +// -------------------------------------------------------------------------- + +#include + +#include +#include + +#include + +#include + +#define EXPAND_OPERATIONS(OPERATION) \ + const unsigned int ls_degree = \ + this->matrix_free.get_dof_handler(parameters.dof_index_ls).get_fe().tensor_degree(); \ + \ + AssertThrow(ls_degree >= 1 && ls_degree <= 4, ExcNotImplemented()); \ + if (ls_degree == 1) \ + OPERATION(1, 0); \ + else if (ls_degree == 2) \ + OPERATION(2, 0); \ + else if (ls_degree == 3) \ + OPERATION(3, 0); \ + else if (ls_degree == 4) \ + OPERATION(4, 0); + +using namespace dealii; + +template +template +void +LevelSetOKZSolverReinitialization::local_reinitialize( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const +{ + const unsigned int concentration_subdivisions = + this->matrix_free.get_dof_handler(parameters.dof_index_ls).get_fe().tensor_degree(); + + const double dtau_inv = std::max(0.95 / (1. / (dim * dim) * this->minimal_edge_length / + concentration_subdivisions), + 1. / (5. * this->time_stepping.step_size())); + + // The second input argument below refers to which constrains should be used, + // 2 means constraints (for LS-function) + FEEvaluation phi(data, + parameters.dof_index_ls, + parameters.quad_index); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(true, true, false); + + VectorizedArray cell_diameter = this->cell_diameters[cell]; + VectorizedArray diffusion = + std::max(make_vectorized_array(this->epsilon_used), + cell_diameter / static_cast(ls_degree)); + + const Tensor<1, dim, VectorizedArray> *normal = + &evaluated_normal[cell * phi.n_q_points]; + for (unsigned int q = 0; q < phi.n_q_points; ++q) + if (!diffuse_only) + { + phi.submit_value(dtau_inv * phi.get_value(q), q); + phi.submit_gradient((diffusion * (normal[q] * phi.get_gradient(q))) * + normal[q], + q); + } + else + { + phi.submit_value(dtau_inv * phi.get_value(q), q); + phi.submit_gradient(phi.get_gradient(q) * diffusion, q); + } + + phi.integrate(true, true); + phi.distribute_local_to_global(dst); + } +} + + + +template +template +void +LevelSetOKZSolverReinitialization::local_reinitialize_rhs( + const MatrixFree &data, + VectorType & dst, + const VectorType &, + const std::pair &cell_range) +{ + // The second input argument below refers to which constrains should be used, + // 2 means constraints (for LS-function) and 4 means constraints_normals + FEEvaluation phi(data, + parameters.dof_index_ls, + parameters.quad_index); + FEEvaluation normals(data, + parameters.dof_index_normal, + parameters.quad_index); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values_plain(this->solution); + phi.evaluate(true, true, false); + + normals.reinit(cell); + normals.read_dof_values_plain(this->normal_vector_field); + normals.evaluate(true, false, false); + + VectorizedArray cell_diameter = this->cell_diameters[cell]; + VectorizedArray diffusion = + std::max(make_vectorized_array(this->epsilon_used), + cell_diameter / static_cast(ls_degree)); + + for (unsigned int q = 0; q < phi.n_q_points; ++q) + if (!diffuse_only) + { + 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()); + evaluated_normal[cell * phi.n_q_points + q] = normal; + } + // take normal as it was for the first reinit step + Tensor<1, dim, VectorizedArray> normal = + evaluated_normal[cell * phi.n_q_points + q]; + phi.submit_gradient(normal * + (0.5 * (1. - phi.get_value(q) * phi.get_value(q)) - + (normal * grad * diffusion)), + q); + } + else + { + phi.submit_gradient(-diffusion * phi.get_gradient(q), q); + } + + phi.integrate(false, true); + phi.distribute_local_to_global(dst); + } +} + + + +template +void +LevelSetOKZSolverReinitialization::reinitialization_vmult( + VectorType & dst, + const VectorType &src, + const bool diffuse_only) const +{ + dst = 0.; + if (diffuse_only) + { +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverReinitialization< \ + dim>::template local_reinitialize, \ + this, \ + dst, \ + src) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + } + else + { +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverReinitialization< \ + dim>::template local_reinitialize, \ + this, \ + dst, \ + src) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + } + + for (unsigned int i = 0; + i < this->matrix_free.get_constrained_dofs(parameters.dof_index_ls).size(); + ++i) + dst.local_element( + this->matrix_free.get_constrained_dofs(parameters.dof_index_ls)[i]) = + preconditioner.get_vector().local_element( + this->matrix_free.get_constrained_dofs(parameters.dof_index_ls)[i]) * + src.local_element( + this->matrix_free.get_constrained_dofs(parameters.dof_index_ls)[i]); +} + + + +template +struct ReinitializationMatrix +{ + ReinitializationMatrix(const LevelSetOKZSolverReinitialization &problem, + const bool diffuse_only) + : problem(problem) + , diffuse_only(diffuse_only) + {} + + void + vmult(VectorType &dst, const VectorType &src) const + { + problem.reinitialization_vmult(dst, src, diffuse_only); + } + + const LevelSetOKZSolverReinitialization &problem; + const bool diffuse_only; +}; + + + +template +void +LevelSetOKZSolverReinitialization::reinitialize(const double dt, + const unsigned int stab_steps, + const unsigned int diff_steps, + const bool) +{ + this->time_stepping.set_time_step(dt); + + // This function assembles and solves for a given profile using the approach + // described in the paper by Olsson, Kreiss, and Zahedi. + + if (evaluated_normal.size() != + this->matrix_free.n_cell_batches() * + this->matrix_free.get_n_q_points(parameters.quad_index)) + evaluated_normal.resize(this->matrix_free.n_cell_batches() * + this->matrix_free.get_n_q_points(parameters.quad_index)); + + std::cout.precision(3); + + // perform several reinitialization steps until we reach the maximum number + // of steps. + // + // TODO: make an adaptive choice of the number of iterations + unsigned actual_diff_steps = diff_steps; + if (this->last_concentration_range.first < -1.02 || + this->last_concentration_range.second > 1.02) + actual_diff_steps += 3; + if (!this->parameters.do_iteration) + this->pcout << (this->time_stepping.now() == this->time_stepping.start() ? " " : + " and ") + << "reinitialize ("; + for (unsigned int tau = 0; tau < actual_diff_steps + stab_steps; tau++) + { + first_reinit_step = (tau == actual_diff_steps); + if (first_reinit_step) + normal_operator.compute_normal(true); + + // compute right hand side + VectorType &rhs = this->system_rhs; + VectorType &increment = this->solution_update; + rhs = 0; + + if (tau < actual_diff_steps) + { +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverReinitialization< \ + dim>::template local_reinitialize_rhs, \ + this, \ + rhs, \ + this->solution) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + } + else + { +#define OPERATION(c_degree, u_degree) \ + this->matrix_free.cell_loop(&LevelSetOKZSolverReinitialization< \ + dim>::template local_reinitialize_rhs, \ + this, \ + rhs, \ + this->solution) + + EXPAND_OPERATIONS(OPERATION); +#undef OPERATION + } + + // solve linear system + { + ReinitializationMatrix matrix(*this, tau < actual_diff_steps); + increment = 0; + + // reduce residual by 1e-6. To obtain good interface shapes, it is + // essential that this tolerance is relative to the rhs + // (ReductionControl steered solver, last argument determines the + // solver) + ReductionControl solver_control(2000, 1e-50, 1e-6); + SolverCG cg(solver_control); + cg.solve(matrix, increment, rhs, preconditioner); + this->constraints.distribute(increment); + if (!this->parameters.do_iteration) + { + if (tau < actual_diff_steps) + this->pcout << "d" << solver_control.last_step(); + else + this->pcout << solver_control.last_step(); + } + } + + this->solution += increment; + this->solution.update_ghost_values(); + + // check residual + const double update_norm = increment.l2_norm(); + if (update_norm < 1e-6) + break; + + if (!this->parameters.do_iteration && tau < actual_diff_steps + stab_steps - 1) + this->pcout << " + "; + } + + if (!this->parameters.do_iteration) + this->pcout << ")" << std::endl << std::flush; + + this->time_stepping.next(); +} + +template class LevelSetOKZSolverReinitialization<2>; +template class LevelSetOKZSolverReinitialization<3>; diff --git a/source/navier_stokes.cc b/source/navier_stokes.cc index 097a94ac..97f93cf9 100644 --- a/source/navier_stokes.cc +++ b/source/navier_stokes.cc @@ -61,6 +61,11 @@ NavierStokes::NavierStokes( parameters.use_simplex_mesh ? std::shared_ptr>(new MappingFE(Simplex::FE_P(1))) : std::shared_ptr>(new MappingQ(3))) + , user_rhs(2) + , solution(2) + , solution_old(2) + , solution_old_old(2) + , solution_update(2) , time_stepping(parameters) , parameters(parameters) , n_mpi_processes(Utilities::MPI::n_mpi_processes(get_communicator(triangulation_in))) @@ -85,6 +90,8 @@ NavierStokes::NavierStokes( , dof_handler_u(triangulation) , dof_handler_p(triangulation) , navier_stokes_matrix(parameters, solution_old, solution_old_old) + , system_rhs(2) + , const_rhs(2) , preconditioner(parameters, *this, triangulation, constraints_u) , dofs_distributed(false) , system_is_setup(false) diff --git a/source/parameters.cc b/source/parameters.cc index 488c132d..e745bcdd 100644 --- a/source/parameters.cc +++ b/source/parameters.cc @@ -576,13 +576,13 @@ FlowParameters::parse_parameters(ParameterHandler &prm) const std::string schem = prm.get("scheme"); if (schem == std::string("implicit_euler")) - time_step_scheme = TimeStepping::implicit_euler; + time_step_scheme = TimeSteppingParameters::Scheme::implicit_euler; else if (schem == std::string("explicit_euler")) - time_step_scheme = TimeStepping::explicit_euler; + time_step_scheme = TimeSteppingParameters::Scheme::explicit_euler; else if (schem == std::string("crank_nicolson")) - time_step_scheme = TimeStepping::crank_nicolson; + time_step_scheme = TimeSteppingParameters::Scheme::crank_nicolson; else if (schem == std::string("bdf_2")) - time_step_scheme = TimeStepping::bdf_2; + time_step_scheme = TimeSteppingParameters::Scheme::bdf_2; else // parameter handler should make sure that we // never end up here diff --git a/source/phase_field.cc b/source/phase_field.cc index 72e69973..0bcb5815 100644 --- a/source/phase_field.cc +++ b/source/phase_field.cc @@ -118,16 +118,9 @@ PhaseFieldSolver::initialize_data_structures() // conditions on open boundaries Functions::ZeroFunction zero_func(1); std::map *> homogeneous_dirichlet; - for (typename std::set::const_iterator it = - this->boundary->fluid_type_plus.begin(); - it != this->boundary->fluid_type_plus.end(); - ++it) - homogeneous_dirichlet[*it] = &zero_func; - for (typename std::set::const_iterator it = - this->boundary->fluid_type_minus.begin(); - it != this->boundary->fluid_type_minus.end(); - ++it) - homogeneous_dirichlet[*it] = &zero_func; + for (const auto &it : this->boundary->fluid_type) + homogeneous_dirichlet[it.first] = &zero_func; + VectorTools::interpolate_boundary_values(this->dof_handler, homogeneous_dirichlet, this->constraints); diff --git a/source/phase_field_local.cc b/source/phase_field_local.cc index 974fbc9b..6cd89e9a 100644 --- a/source/phase_field_local.cc +++ b/source/phase_field_local.cc @@ -192,7 +192,7 @@ PhaseFieldSolver::local_residual( w_values[q] = c_val; val += this->time_stepping.weight_old() * inv_time_weight * old_c_values.get_value(q); - if (this->time_stepping.scheme() == TimeStepping::bdf_2 && + if (this->time_stepping.scheme() == TimeSteppingParameters::Scheme::bdf_2 && this->time_stepping.step_no() > 1) val += this->time_stepping.weight_old_old() * inv_time_weight * old_old_c_values.get_value(q); diff --git a/source/time_stepping.cc b/source/time_stepping.cc index 8dcca9e9..c6247514 100644 --- a/source/time_stepping.cc +++ b/source/time_stepping.cc @@ -21,7 +21,6 @@ TimeStepping::TimeStepping(const FlowParameters ¶meters) : start_val(parameters.start_time) , final_val(parameters.end_time) - , tolerance_val(parameters.time_step_tolerance) , scheme_val(parameters.time_step_scheme) , start_step_val(parameters.time_step_size_start) , max_step_val(parameters.time_step_size_max) @@ -40,26 +39,64 @@ TimeStepping::TimeStepping(const FlowParameters ¶meters) { now_val = start_val; prev_val = start_val; - if (scheme_val == implicit_euler) + if (scheme_val == TimeSteppingParameters::Scheme::implicit_euler) { tau1_val = 1.; tau2_val = 0.; } - else if (scheme_val == explicit_euler) + else if (scheme_val == TimeSteppingParameters::Scheme::explicit_euler) { tau1_val = 0.; tau2_val = 1.; } - else if (scheme_val == crank_nicolson) + else if (scheme_val == TimeSteppingParameters::Scheme::crank_nicolson) tau1_val = tau2_val = .5; - else if (scheme_val == bdf_2) + else if (scheme_val == TimeSteppingParameters::Scheme::bdf_2) { tau1_val = 1.; tau2_val = 0.; } } - +TimeStepping::TimeStepping(const TimeSteppingParameters ¶meters) + : start_val(parameters.start_time) + , final_val(parameters.end_time) + , scheme_val(parameters.time_step_scheme) + , start_step_val(parameters.time_step_size_start) + , max_step_val(parameters.time_step_size_max) + , min_step_val(parameters.time_step_size_min) + , current_step_val(start_step_val) + , last_step_val(0.) + , step_val(start_step_val) + , weight_val(1. / start_step_val) + , weight_old_val(-1.) + , weight_old_old_val(0.) + , factor_extrapol_old(0.) + , factor_extrapol_old_old(0.) + , step_no_val(0) + , at_end_val(false) + , weight_changed(true) +{ + now_val = start_val; + prev_val = start_val; + if (scheme_val == TimeSteppingParameters::Scheme::implicit_euler) + { + tau1_val = 1.; + tau2_val = 0.; + } + else if (scheme_val == TimeSteppingParameters::Scheme::explicit_euler) + { + tau1_val = 0.; + tau2_val = 1.; + } + else if (scheme_val == TimeSteppingParameters::Scheme::crank_nicolson) + tau1_val = tau2_val = .5; + else if (scheme_val == TimeSteppingParameters::Scheme::bdf_2) + { + tau1_val = 1.; + tau2_val = 0.; + } +} void TimeStepping::restart() @@ -91,7 +128,7 @@ TimeStepping::next() if (now_val != start()) { last_step_val = current_step_val; - if (scheme_val == bdf_2 && step_no_val == 1) + if (scheme_val == TimeSteppingParameters::Scheme::bdf_2 && step_no_val == 1) s = step_val; if (s > max_step_val) @@ -117,7 +154,7 @@ TimeStepping::next() { double new_weight; - if (scheme_val == bdf_2 && now_val != start()) + if (scheme_val == TimeSteppingParameters::Scheme::bdf_2 && now_val != start()) { new_weight = ((2. * current_step_val + last_step_val) / (current_step_val * (current_step_val + last_step_val))); @@ -165,20 +202,19 @@ std::string TimeStepping::name() const { std::string result; - if (scheme_val == implicit_euler) + if (scheme_val == TimeSteppingParameters::Scheme::implicit_euler) result = std::string("ImplEuler"); - else if (scheme_val == explicit_euler) + else if (scheme_val == TimeSteppingParameters::Scheme::explicit_euler) result = std::string("ExplEuler"); - else if (scheme_val == crank_nicolson) + else if (scheme_val == TimeSteppingParameters::Scheme::crank_nicolson) result = std::string("CrankNicolson"); - else if (scheme_val == bdf_2) + else if (scheme_val == TimeSteppingParameters::Scheme::bdf_2) result = std::string("BDF-2"); return result; } - -TimeStepping::Scheme +TimeSteppingParameters::Scheme TimeStepping::scheme() const { return scheme_val; diff --git a/source/two_phase_base.cc b/source/two_phase_base.cc index 1d6de804..b1422e89 100644 --- a/source/two_phase_base.cc +++ b/source/two_phase_base.cc @@ -65,7 +65,11 @@ TwoPhaseBaseAlgorithm::TwoPhaseBaseAlgorithm( const std::shared_ptr> fe_in, Triangulation & tria_in, TimerOutput * timer_in) - : pcout(std::cout, Utilities::MPI::this_mpi_process(get_communicator(tria_in)) == 0) + : solution_update(2) + , solution(2) + , solution_old(2) + , solution_old_old(2) + , pcout(std::cout, Utilities::MPI::this_mpi_process(get_communicator(tria_in)) == 0) , timer((timer_in == 0 ? new TimerOutput(pcout, parameters_in.output_wall_times ? TimerOutput::summary : @@ -77,6 +81,7 @@ TwoPhaseBaseAlgorithm::TwoPhaseBaseAlgorithm( , navier_stokes(parameters_in, triangulation, timer.get(), this->boundary) , fe(fe_in) , dof_handler(triangulation) + , system_rhs(2) , time_stepping(navier_stokes.time_stepping) , parameters(navier_stokes.get_parameters()) , epsilon_used(0)