From 171e600d99e2c00c9f384ff47cebcc6bb4d3ca72 Mon Sep 17 00:00:00 2001 From: landinjm Date: Mon, 27 Jan 2025 14:28:59 -0500 Subject: [PATCH] nonuniform dirichlet --- applications/multigrid_test/CMakeLists.txt | 2 +- applications/multigrid_test/ICs_and_BCs.cc | 24 ++++ applications/multigrid_test/equations.cc | 6 +- applications/multigrid_test/parameters.prm | 4 +- .../boundary_conditions/constraint_handler.h | 22 ++-- .../nonuniform_dirichlet.h | 115 ++++++++++++++++++ include/core/pde_problem.h | 14 ++- include/core/triangulation_handler.h | 69 +++++++++++ src/core/variable_container.cc | 4 +- tests/core/field_solve_type.cc | 2 +- 10 files changed, 241 insertions(+), 21 deletions(-) create mode 100644 applications/multigrid_test/ICs_and_BCs.cc create mode 100644 include/core/boundary_conditions/nonuniform_dirichlet.h diff --git a/applications/multigrid_test/CMakeLists.txt b/applications/multigrid_test/CMakeLists.txt index 19e5db0f9..7174ecfd9 100644 --- a/applications/multigrid_test/CMakeLists.txt +++ b/applications/multigrid_test/CMakeLists.txt @@ -15,7 +15,7 @@ include_directories(${CMAKE_SOURCE_DIR}/../../src) include_directories(${CMAKE_SOURCE_DIR}) # Set the location of the main.cc file -set(TARGET_SRC "${CMAKE_SOURCE_DIR}/main.cc" "${CMAKE_SOURCE_DIR}/equations.cc") +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/main.cc" "${CMAKE_SOURCE_DIR}/equations.cc" "${CMAKE_SOURCE_DIR}/ICs_and_BCs.cc") # Check if there has been updates to main library set(dir ${PROJECT_SOURCE_DIR}/../..) diff --git a/applications/multigrid_test/ICs_and_BCs.cc b/applications/multigrid_test/ICs_and_BCs.cc new file mode 100644 index 000000000..cac330149 --- /dev/null +++ b/applications/multigrid_test/ICs_and_BCs.cc @@ -0,0 +1,24 @@ +#include +#include +#include + +template +void +customNonuniformDirichlet::set_nonuniform_dirichlet( + [[maybe_unused]] const uint &index, + [[maybe_unused]] const uint &boundary_id, + [[maybe_unused]] const uint &component, + [[maybe_unused]] const dealii::Point &point, + [[maybe_unused]] double &scalar_value, + [[maybe_unused]] double &vector_component_value) const +{ + if (index == 0) + { + if (boundary_id == 0 || boundary_id == 1) + { + scalar_value = std::sin(point[1] * M_PI); + } + } +} + +INSTANTIATE_UNI_TEMPLATE(customNonuniformDirichlet) \ No newline at end of file diff --git a/applications/multigrid_test/equations.cc b/applications/multigrid_test/equations.cc index b2b8cecb9..0e708d666 100644 --- a/applications/multigrid_test/equations.cc +++ b/applications/multigrid_test/equations.cc @@ -11,7 +11,7 @@ customAttributeLoader::loadVariableAttributes() set_variable_type(0, SCALAR); set_variable_equation_type(0, TIME_INDEPENDENT); - set_dependencies_value_term_RHS(0, "grad(phi)"); + set_dependencies_value_term_RHS(0, ""); set_dependencies_gradient_term_RHS(0, "grad(phi)"); set_dependencies_value_term_LHS(0, ""); set_dependencies_gradient_term_LHS(0, "grad(change(phi))"); @@ -36,10 +36,6 @@ customPDE::compute_nonexplicit_RHS( { scalarGrad phix = variable_list.get_scalar_gradient(0); - scalarValue coefficient = 1.0 / (0.05 + 2.0 * q_point_loc.square()); - - variable_list.set_scalar_value_term(0, - make_vectorized_array(1.0) / coefficient); variable_list.set_scalar_gradient_term(0, -phix); } diff --git a/applications/multigrid_test/parameters.prm b/applications/multigrid_test/parameters.prm index d0250356e..61bcecf66 100644 --- a/applications/multigrid_test/parameters.prm +++ b/applications/multigrid_test/parameters.prm @@ -8,9 +8,9 @@ set Subdivisions X = 1 set Subdivisions Y = 1 set Subdivisions Z = 1 set Refine factor = 5 -set Element degree = 6 +set Element degree = 1 -set Boundary condition for variable phi = DIRICHLET: 0.0 +set Boundary condition for variable phi = NON_UNIFORM_DIRICHLET subsection Pinning point: phi set value = 2147483647 diff --git a/include/core/boundary_conditions/constraint_handler.h b/include/core/boundary_conditions/constraint_handler.h index 57bc63d7f..5c7ba3dbf 100644 --- a/include/core/boundary_conditions/constraint_handler.h +++ b/include/core/boundary_conditions/constraint_handler.h @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -24,8 +25,8 @@ class constraintHandler * \brief Constructor. */ constraintHandler(const boundaryParameters &_boundary_parameters, - const uint &global_variable_index, - const fieldType &field_type); + const uint &_global_variable_index, + const fieldType &_field_type); /** * \brief Getter function for constraints (constant reference). @@ -55,11 +56,11 @@ class constraintHandler std::unique_ptr> constraints; - const boundaryParameters &boundary_parameters; + boundaryParameters boundary_parameters; - const uint &global_variable_index; + uint global_variable_index; - const fieldType &field_type; + fieldType field_type; }; template @@ -104,6 +105,7 @@ constraintHandler::make_constraints(const dealii::MappingQ1 &mapping, { if (boundary_type == boundaryType::NATURAL) { + // Do nothing because they are naturally enforced. continue; } else if (boundary_type == boundaryType::DIRICHLET) @@ -119,6 +121,7 @@ constraintHandler::make_constraints(const dealii::MappingQ1 &mapping, } else if (boundary_type == boundaryType::PERIODIC) { + // Do nothing because this is enforced in triangulation_handler. continue; } else if (boundary_type == boundaryType::NEUMANN) @@ -127,8 +130,13 @@ constraintHandler::make_constraints(const dealii::MappingQ1 &mapping, } else if (boundary_type == boundaryType::NON_UNIFORM_DIRICHLET) { - Assert(false, - FeatureNotImplemented("Nonuniform dirichlet boundary conditions")); + dealii::VectorTools::interpolate_boundary_values( + mapping, + dof_handler, + boundary_id, + nonuniformDirichlet(global_variable_index, + boundary_id), + *constraints); continue; } else if (boundary_type == boundaryType::NON_UNIFORM_NEUMANN) diff --git a/include/core/boundary_conditions/nonuniform_dirichlet.h b/include/core/boundary_conditions/nonuniform_dirichlet.h new file mode 100644 index 000000000..eeb0d44be --- /dev/null +++ b/include/core/boundary_conditions/nonuniform_dirichlet.h @@ -0,0 +1,115 @@ +#ifndef nonuniform_dirichlet_h +#define nonuniform_dirichlet_h + +#include +#include + +#include + +/** + * \brief Forward declaration of user-facing implementation + */ +template +class customNonuniformDirichlet; + +/** + * \brief Function for user-implemented nonuniform dirichlet boundary condition. + */ +template +class nonuniformDirichlet : public dealii::Function +{ +public: + /** + * \brief Constructor. + */ + nonuniformDirichlet(const uint &_index, const uint &_boundary_id) + : dealii::Function((field_type == fieldType::VECTOR) ? dim : 1) + , index(_index) + , boundary_id(_boundary_id) {}; + + /** + * \brief Scalar value. + */ + double + value(const dealii::Point &p, const unsigned int component = 0) const override; + + /** + * \brief Vector value. + */ + void + vector_value(const dealii::Point &p, dealii::Vector &value) const override; + +private: + const uint index; + + const uint boundary_id; + + customNonuniformDirichlet custom_nonuniform_dirichlet; +}; + +template +inline double +nonuniformDirichlet::value( + const dealii::Point &p, + [[maybe_unused]] const unsigned int component) const +{ + // Initialize passed variables to zero + double scalar_value = 0.0; + dealii::Vector vector_value(dim); + + // Pass variables to user-facing function to evaluate + custom_nonuniform_dirichlet + .set_nonuniform_dirichlet(index, boundary_id, 0, p, scalar_value, vector_value(0)); + + return scalar_value; +} + +template +inline void +nonuniformDirichlet::vector_value(const dealii::Point &p, + dealii::Vector &value) const +{ + // Initialize passed variables to zero + double scalar_value = 0.0; + dealii::Vector vector_value(dim); + + // Pass variables to user-facing function to evaluate + for (uint i = 0; i < dim; i++) + { + custom_nonuniform_dirichlet.set_nonuniform_dirichlet(index, + boundary_id, + i, + p, + scalar_value, + vector_value(i)); + } + + value = vector_value; +} + +/** + * \brief User-facing implementation of nonuniform boundary conditions + */ +template +class customNonuniformDirichlet +{ +public: + /** + * \brief Constructor. + */ + customNonuniformDirichlet() = default; + + /** + * \brief Function that passes the value/vector and point that are set in the nonuniform + * dirichlet. + */ + void + set_nonuniform_dirichlet(const uint &index, + const uint &boundary_id, + const uint &component, + const dealii::Point &point, + double &scalar_value, + double &vector_component_value) const; +}; + +#endif \ No newline at end of file diff --git a/include/core/pde_problem.h b/include/core/pde_problem.h index 596cb17df..0b481fbbe 100644 --- a/include/core/pde_problem.h +++ b/include/core/pde_problem.h @@ -176,18 +176,26 @@ PDEProblem::init_system() system_matrix.initialize_dof_vector(solution); system_matrix.initialize_dof_vector(residual); + // Apply boundary conditions to the solution vector as an initial guess + constraint_handler.get_constraints().distribute(solution); + const uint nlevels = triangulation_handler.get_n_global_levels(); mg_matrices.resize(0, nlevels - 1, user_inputs); - const std::set dirichlet_boundary_ids = {0}; + // This should correspond to the userinput dirichlets + const std::set dirichlet_boundary_ids = {0, 1, 2, 3}; + mg_constrained_dofs.initialize(dof_handler); mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary_ids); for (uint level = 0; level < nlevels; ++level) { + // Setup the constraint set for the multigrid level dealii::AffineConstraints level_constraints( dof_handler.locally_owned_mg_dofs(level), dealii::DoFTools::extract_locally_relevant_level_dofs(dof_handler, level)); + + // Loop through boundary indices and apply constraints for (const dealii::types::global_dof_index dof_index : mg_constrained_dofs.get_boundary_indices(level)) { @@ -282,7 +290,7 @@ PDEProblem::solve_increment() system_matrix.compute_residual(residual, solution); - dealii::SolverControl solver_control(100, 1e-12 * residual.l2_norm()); + dealii::SolverControl solver_control(1000, 1e-12 * residual.l2_norm()); dealii::SolverCG> cg(solver_control); Timer time; @@ -291,7 +299,7 @@ PDEProblem::solve_increment() try { change_in_solution = 0.0; - cg.solve(system_matrix, change_in_solution, residual, preconditioner); + cg.solve(system_matrix, change_in_solution, residual, dealii::IdentityMatrix()); constraint_handler.get_constraints().set_zero(change_in_solution); } catch (...) diff --git a/include/core/triangulation_handler.h b/include/core/triangulation_handler.h index 512d7a554..5c0dbb8fb 100644 --- a/include/core/triangulation_handler.h +++ b/include/core/triangulation_handler.h @@ -2,10 +2,13 @@ #define triangulation_handler_h #include +#include #include +#include #include #include +#include /** * \brief This class handlers the generation and manipulation of triangulations. @@ -42,7 +45,20 @@ class triangulationHandler void generate_mesh(const userInputParameters &user_inputs); + /** + * \brief Export triangulation to vtk. + */ + void + export_triangulation_as_vtk(const std::string &filename) const; + private: + /** + * \brief Mark the domain ids on the triangulation to get the proper mapping of + * specified boundary conditions. + */ + void + mark_boundaries(const userInputParameters &user_inputs) const; + std::unique_ptr triangulation; }; @@ -88,8 +104,61 @@ triangulationHandler::generate_mesh(const userInputParameters &user_in dealii::Point(), dealii::Point(user_inputs.spatial_discretization.domain_size)); + // Mark boundaries. This is done before global refinement to reduce the number of cells + // we have to loop through. + mark_boundaries(user_inputs); + + // Output triangulation to vtk if in debug mode +#ifdef DEBUG + export_triangulation_as_vtk("triangulation"); +#endif + // Global refinement triangulation->refine_global(user_inputs.spatial_discretization.refine_factor); } +template +void +triangulationHandler::export_triangulation_as_vtk(const std::string &filename) const +{ + dealii::GridOut grid_out; + std::ofstream out(filename + ".vtk"); + grid_out.write_vtk(*triangulation, out); + std::cout << "Triangulation written to " << filename << ".vtk\n"; +} + +template +void +triangulationHandler::mark_boundaries( + const userInputParameters &user_inputs) const +{ + double tolerance = 1e-12; + + // Loop through the cells + for (const auto &cell : triangulation->active_cell_iterators()) + { + // Mark the faces (faces_per_cell = 2*dim) + for (uint face_number = 0; face_number < dealii::GeometryInfo::faces_per_cell; + ++face_number) + { + // Direction for quad and hex cells + uint direction = std::floor(face_number / 2); + + // Mark the boundary id for x=0, y=0, z=0 + if (std::fabs(cell->face(face_number)->center()(direction) - 0) < tolerance) + { + cell->face(face_number)->set_boundary_id(face_number); + } + // Mark the boundary id for x=max, y=max, z=max + else if (std::fabs( + cell->face(face_number)->center()(direction) - + (user_inputs.spatial_discretization.domain_size[direction])) < + tolerance) + { + cell->face(face_number)->set_boundary_id(face_number); + } + } + } +} + #endif \ No newline at end of file diff --git a/src/core/variable_container.cc b/src/core/variable_container.cc index 36bcc6566..0dbe9dbf0 100644 --- a/src/core/variable_container.cc +++ b/src/core/variable_container.cc @@ -271,7 +271,7 @@ variableContainer::reinit_and_eval( auto *scalar_FEEval_ptr = scalar_vars_map.at(dependency_index).at(dependency_type).get(); scalar_FEEval_ptr->reinit(cell); - scalar_FEEval_ptr->read_dof_values(src); + scalar_FEEval_ptr->read_dof_values_plain(src); scalar_FEEval_ptr->evaluate(flags); } else @@ -279,7 +279,7 @@ variableContainer::reinit_and_eval( auto *vector_FEEval_ptr = vector_vars_map.at(dependency_index).at(dependency_type).get(); vector_FEEval_ptr->reinit(cell); - vector_FEEval_ptr->read_dof_values(src); + vector_FEEval_ptr->read_dof_values_plain(src); vector_FEEval_ptr->evaluate(flags); } } diff --git a/tests/core/field_solve_type.cc b/tests/core/field_solve_type.cc index 9dd859c6d..e7ae14371 100644 --- a/tests/core/field_solve_type.cc +++ b/tests/core/field_solve_type.cc @@ -80,7 +80,7 @@ TEST_CASE("Field solve types") * weak form this becomes * ∇w⋅(∇u/√(1+|∇u|^2))+∇w⋅(δ∇u/√(1+|∇u|^2))-∇w⋅((∇u⋅δ∇u)∇u/√(1+|∇u|^2))=0. * - * The third equuation is a simple Poisson equation Δu=0. + * The third equation is a simple Poisson equation Δu=0. * * The fourth equation is the steady-state version of the Allen-Cahn equation * u^3-u-γΔu=0. In the weak form with γ=1 this becomes w(u^3+δu^3)+w(u-δu)+∇w(∇u+δ∇u)=0.