diff --git a/applications/allen_cahn_conserved/CMakeLists.txt b/applications/allen_cahn_conserved/CMakeLists.txt new file mode 100644 index 000000000..c1ecb3245 --- /dev/null +++ b/applications/allen_cahn_conserved/CMakeLists.txt @@ -0,0 +1,73 @@ +## +# CMake script for the PRISMS-PF applications +## + +cmake_minimum_required(VERSION 3.28) +cmake_policy(VERSION 3.28) + +# Name the application +set(APPLICATION_NAME "allen_cahn_conserved") + +# Create a project for the application +project(${APPLICATION_NAME} CXX) + +# Set the build type to debug is not specified +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Build type" FORCE) +endif() + +# Export compile commands for IDEs +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +# Find the PRISMS-PF package +find_package( + PRISMS-PF + REQUIRED + HINTS + ${PRISMS_PF_DIR} + $ENV{PRISMS_PF_DIR} +) + +# Declare the application executable +add_executable(${APPLICATION_NAME}) + +# Rename the output executable to main +set_target_properties( + ${APPLICATION_NAME} + PROPERTIES + OUTPUT_NAME + main +) + +# Add sources +target_sources( + ${APPLICATION_NAME} + PRIVATE + main.cc + custom_pde.h +) + +# Require C++20 +# TODO: Remove this +target_compile_features(${APPLICATION_NAME} PRIVATE cxx_std_20) +set_target_properties( + ${APPLICATION_NAME} + PROPERTIES + CXX_STANDARD + 20 + CXX_STANDARD_REQUIRED + ON + CXX_EXTENSIONS + OFF +) + +# Link PRISMS-PF +# If we're in Debug configuration and PRISMS_PF_BUILD_TYPE is DebugRelease, link to the prisms_pf_debug target, +# which is only available in prisms_pf DebugRelease configuration. +# Otherwise link to the prisms_pf target, which may be debug or release +target_link_libraries( + ${APPLICATION_NAME} + PRIVATE + $,$>,PRISMS-PF::prisms_pf_debug, + PRISMS-PF::prisms_pf> +) diff --git a/applications/allen_cahn_conserved/allen_cahn_conserved.md b/applications/allen_cahn_conserved/allen_cahn_conserved.md new file mode 100644 index 000000000..a7d072655 --- /dev/null +++ b/applications/allen_cahn_conserved/allen_cahn_conserved.md @@ -0,0 +1,118 @@ +# PRISMS PhaseField: Globally Conserved Allen-Cahn Dynamics + +This application performs a phase field simulation of Allen-Cahn dynamics subject to global (as opposed to local) conservation. Global conservation implies that $\int_\Omega c \,\mathrm{d}V$ is constant in time. + +Consider a free energy expression of the form: + +$$ +\begin{equation} +F[c] = \int_{\Omega} f(c) + \frac{\kappa}{2} \nabla c \cdot \nabla c \,\mathrm{d} V, +\end{equation} +$$ + +where $f(c)$ is the chemical free energy density, $\kappa$ is the gradient coefficient, and $c$ is the structural order parameter, which subject to a global constraint + +$$ +\begin{equation} +\int_{\Omega} c \,\mathrm{d} V = \mathrm{constant}. +\end{equation} +$$ + +The kinetics is given by the Allen-Cahn equation, which is the gradient flow of the functional $F[c]$: + +$$ +\begin{equation} +\frac{\partial c}{\partial t} = - M \mu, +\end{equation} +$$ + +where $M$ is the mobility and $\mu$ is the chemical potential defined as + +$$ +\begin{equation} +\mu = \frac{\delta F}{\delta c} = f_{,c}(c) - \kappa \nabla^2 c. +\end{equation} +$$ + +To conserve $c$, we can introduce a global Lagrange multiplier $\lambda(t)$ + +$$ +\begin{equation} +\frac{\partial c}{\partial t} = - M (\mu - \lambda(t)). +\end{equation} +$$ + +The Lagrange multiplier can be solved from the conservation of $c$ (Eq. 2) as + +$$ +\begin{equation} +\lambda(t) = \frac{1}{V} \int_{\Omega} \mu \,\mathrm{d} V, \quad V = \int_{\Omega} \,\mathrm{d} V. +\end{equation} +$$ + +The governing equations are given by Eqs. 5, 4, and 6. + + +## Time discretization +Considering forward Euler explicit time stepping, we have the time discretized kinetics equations: + +$$ +\begin{equation} +c^{n+1} = c^{n} - \Delta t M(\mu^n-\lambda^n) +\end{equation} +$$ + +and + +$$ +\begin{equation} +\mu^{n} = f_{,c}^n - \kappa \nabla^2 c^n. +\end{equation} +$$ + +## Weak formulation +In the weak formulation, considering an arbitrary variation $w$, the above equation can be expressed as a residual equations: + +$$ +\begin{equation} +\int_{\Omega} w c^{n+1} \,\mathrm{d}V = \int_{\Omega} w \left( c^{n} - \Delta t M(\mu^n-\lambda^n) \right) \,\mathrm{d}V +\end{equation} +$$ + +$$ +\begin{equation} +r_{c} = c^{n} - \Delta t M(\mu^n-\lambda^n) +\end{equation} +$$ + +and + +$$ +\begin{align} +\int_{\Omega} w \mu^{n} \,\mathrm{d}V &= \int_{\Omega} w f_{,c}^n - w \kappa \nabla^2 c^{n} \,\mathrm{d}V \\ +&= \int_{\Omega} w ~ (f_{,c}^{n}) + \nabla w \cdot (\kappa \nabla c^{n}) \,\mathrm{d}V, +\end{align} +$$ + +$$ +\begin{equation} +r_{\mu} = f_{,\eta}^{n} +\end{equation} +$$ + +$$ +\begin{equation} +r_{\mu x} = \kappa \nabla \eta^{n} +\end{equation} +$$ + +where the Lagrange multiplier $\lambda^n$ for time step $n$ is calculated as + +$$ +\begin{align} +\lambda^n& = \frac{1}{V}\int_\Omega \mu^n \,\mathrm{d}V. +\end{align} +$$ + +The above values of $r_{c}$, $r_{\mu}$, and $r_{\mu x}$ are used in `compute_rhs()` in the following source file: +`applications/allen_cahn_conserved/custom_pde.cc` diff --git a/applications/allen_cahn_conserved/custom_pde.h b/applications/allen_cahn_conserved/custom_pde.h new file mode 100644 index 000000000..5c6f102fa --- /dev/null +++ b/applications/allen_cahn_conserved/custom_pde.h @@ -0,0 +1,167 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include + +#include + +#include + +#include + +PRISMS_PF_BEGIN_NAMESPACE + +template +class CustomPDE : public PDEOperatorBase +{ +public: + using ScalarValue = dealii::VectorizedArray; + using ScalarGrad = dealii::Tensor<1, dim, ScalarValue>; + using ScalarHess = dealii::Tensor<2, dim, ScalarValue>; + using VectorValue = dealii::Tensor<1, dim, ScalarValue>; + using VectorGrad = dealii::Tensor<2, dim, ScalarValue>; + using VectorHess = dealii::Tensor<3, dim, ScalarValue>; + using PDEOperatorBase::get_user_inputs; + using PDEOperatorBase::get_pf_tools; + + /** + * @brief Constructor. + */ + explicit CustomPDE(const UserInputParameters &_user_inputs, + PhaseFieldTools &_pf_tools) + : PDEOperatorBase(_user_inputs, _pf_tools) + , mobility(get_user_inputs().user_constants.get_double("mobility")) + , kappa(get_user_inputs().user_constants.get_double("kappa")) + {} + + void + set_lambda(number v) + { + lambda = v; + } + +private: + number mobility; + number kappa; + number lambda; + + void + set_initial_condition([[maybe_unused]] const unsigned int &index, + [[maybe_unused]] const unsigned int &component, + [[maybe_unused]] const dealii::Point &point, + [[maybe_unused]] number &scalar_value, + [[maybe_unused]] number &vector_component_value) const override + { + const dealii::Tensor<1, dim> &mesh_size = + get_user_inputs().spatial_discretization.rectangular_mesh.size; + constexpr number center[12][3] = { + {0.1, 0.3, 0}, + {0.8, 0.7, 0}, + {0.5, 0.2, 0}, + {0.4, 0.4, 0}, + {0.3, 0.9, 0}, + {0.8, 0.1, 0}, + {0.9, 0.5, 0}, + {0.0, 0.1, 0}, + {0.1, 0.6, 0}, + {0.5, 0.6, 0}, + {1, 1, 0}, + {0.7, 0.95, 0} + }; + constexpr number rad[12] = {12, 14, 19, 16, 11, 12, 17, 15, 20, 10, 11, 14}; + + number dist = 0.0; + number sdf = std::numeric_limits::max(); + for (unsigned int i = 0; i < 12; i++) + { + dist = 0.0; + for (unsigned int dir = 0; dir < dim; dir++) + { + const number comp_diff = point[dir] - center[i][dir] * mesh_size[dir]; + dist += comp_diff * comp_diff; + } + dist = std::sqrt(dist) - rad[i]; + sdf = std::min(sdf, dist); + } + scalar_value = 0.5 * (1.0 - std::tanh(sdf / 1.5)); + } + + void + compute_rhs(FieldContainer &variable_list, + const SimulationTimer &sim_timer, + unsigned int solve_block_id) const override + { + if (solve_block_id == 0) // explicit c + { + ScalarValue c = variable_list.template get_value(0); + ScalarValue mu = variable_list.template get_value(1); + + ScalarValue eq_c = c - sim_timer.get_timestep() * mobility * (mu - lambda); + + variable_list.set_value_term(0, eq_c); + } + else if (solve_block_id == 1) // mu + { + ScalarValue c = variable_list.template get_value(0); + ScalarGrad c_grad = variable_list.template get_gradient(0); + ScalarValue df_well = 4.0 * c * (c - 1.0) * (c - 0.5); + + variable_list.set_value_term(1, df_well); + variable_list.set_gradient_term(1, kappa * c_grad); + } + else if (solve_block_id == 2) // custom + { + } + else if (solve_block_id == 3) // postprocess + { + ScalarValue c = variable_list.template get_value(0); + ScalarGrad c_grad = variable_list.template get_gradient(0); + ScalarValue f_chem = c * c * (1.0 - c) * (1.0 - c); + ScalarValue f_grad = 0.5 * kappa * c_grad.norm_square(); + + variable_list.set_value_term(3, c_grad.norm()); + variable_list.set_value_term(4, f_chem + f_grad); + } + } +}; + +template +class MyCustomSolver : public SolverBase +{ +private: + CustomPDE &custom_pde; + +public: + MyCustomSolver(const SolveBlock &solve_block, + const SolveContext &solve_context, + CustomPDE &_custom_pde) + : SolverBase(solve_block, solve_context) + , custom_pde(_custom_pde) + {} + + void + init(const std::list &all_dependency_sets) override + { + SolverBase::init(all_dependency_sets); + } + + void + solve_impl() override + { + number mu_int = Integrator::template integrate<0>( + this->solve_context->get_dof_manager().get_field_dof_handler(1), + this->solve_context->get_solution_indexer().get_solution_vector(1))[0]; + + const dealii::Tensor<1, dim> &mesh_size = + this->solve_context->get_user_inputs().spatial_discretization.rectangular_mesh.size; + number V = mesh_size[0] * mesh_size[1]; + + custom_pde.set_lambda(mu_int / V); + } + + void + update() override + {} +}; + +PRISMS_PF_END_NAMESPACE diff --git a/applications/allen_cahn_conserved/main.cc b/applications/allen_cahn_conserved/main.cc new file mode 100644 index 000000000..c676b381a --- /dev/null +++ b/applications/allen_cahn_conserved/main.cc @@ -0,0 +1,88 @@ +// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan +// SPDX-License-Identifier: GNU Lesser General Public Version 2.1 + +#include "custom_pde.h" + +#include +#include +#include +#include + +using namespace prismspf; + +int +main(int argc, char *argv[]) +{ + // Initialize MPI + prismspf::MPIInitFinalize mpi_init(argc, argv); + + // Parse the command line options (if there are any) to get the name of the input + // file + ParseCMDOptions cli_options(argc, argv); + + constexpr unsigned int dim = 2; + constexpr unsigned int degree = 1; + + std::vector field_attributes = { + FieldAttributes("c", Scalar), + FieldAttributes("mu", Scalar), + FieldAttributes("dummy", Scalar), + FieldAttributes("mg_n", Scalar), + FieldAttributes("f_tot", Scalar), + }; + + SolveBlock c_block; + c_block.id = 0; + c_block.solve_type = Explicit; + c_block.solve_timing = Primary; + c_block.field_indices = {0}; + c_block.dependencies_rhs = + make_dependency_set(field_attributes, {"old_1(c)", "old_1(mu)"}); + + SolveBlock mu_block; + mu_block.id = 1; + mu_block.solve_type = Explicit; + mu_block.solve_timing = Uninitialized; + mu_block.field_indices = {1}; + mu_block.dependencies_rhs = make_dependency_set(field_attributes, {"c", "grad(c)"}); + + SolveBlock cs_block; + cs_block.id = 2; + cs_block.solve_type = CustomSolver; + cs_block.solve_timing = Uninitialized; + cs_block.field_indices = {2}; + cs_block.dependencies_rhs = make_dependency_set(field_attributes, {"mu"}); + + SolveBlock pp_block; + pp_block.id = 3; + pp_block.solve_type = Explicit; + pp_block.solve_timing = PostProcess; + pp_block.field_indices = {3, 4}; + pp_block.dependencies_rhs = make_dependency_set(field_attributes, {"c", "grad(c)"}); + + UserInputParameters user_inputs(cli_options.get_parameters_filename()); + PhaseFieldTools pf_tools; + CustomPDE pde_operator(user_inputs, pf_tools); + + // setup CustomSolver + using CustomFactory = std::function>( + const SolveBlock &, + const SolveContext &)>; + + cs_block.custom_solver_factory = CustomFactory( + [&pde_operator](const SolveBlock &b, const SolveContext &c) + { + return std::make_shared>(b, c, pde_operator); + }); + + std::vector solve_blocks({c_block, mu_block, cs_block, pp_block}); + + Problem problem(field_attributes, + solve_blocks, + user_inputs, + pf_tools, + pde_operator); + problem.solve(); + + return 0; +} diff --git a/applications/allen_cahn_conserved/parameters.prm b/applications/allen_cahn_conserved/parameters.prm new file mode 100644 index 000000000..2860b7328 --- /dev/null +++ b/applications/allen_cahn_conserved/parameters.prm @@ -0,0 +1,29 @@ +subsection Rectangular mesh + set x size = 100.0 + set y size = 100.0 + set z size = 100.0 +end + +set global refinement = 8 +set max refinement = 8 +set min refinement = 3 +set mesh adaptivity = true +set remeshing period = 800 + +subsection refinement criterion: 0 + set type = value + set value lower bound = 0.01 + set value upper bound = 0.99 + set variables = c +end + +set time step = 1.0e-2 +set number steps = 10000 + +subsection output + set condition = EQUAL_SPACING + set number = 5 +end + +set Model constant mobility = 1.0, DOUBLE +set Model constant kappa = 2.0, DOUBLE diff --git a/include/prismspf/core/dependency_extents.h b/include/prismspf/core/dependency_extents.h index 6b855a600..e4ad82f05 100644 --- a/include/prismspf/core/dependency_extents.h +++ b/include/prismspf/core/dependency_extents.h @@ -35,9 +35,9 @@ struct DependencyExtents // TODO: update this to account multigrid levels. Will need more // than std::list DependencyExtents(const std::set &field_indices, - const std::list &all_dependeny_sets) + const std::list &all_dependency_sets) { - for (const DependencyMap &dependency_set : all_dependeny_sets) + for (const DependencyMap &dependency_set : all_dependency_sets) { for (unsigned int field_index : field_indices) { diff --git a/include/prismspf/core/solve_block.h b/include/prismspf/core/solve_block.h index 4d8c14698..a55f4930e 100644 --- a/include/prismspf/core/solve_block.h +++ b/include/prismspf/core/solve_block.h @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -136,6 +137,11 @@ class SolveBlock return id < other.id; } + /** + * @brief Type-erased factory function for instantiating custom solvers. + */ + std::any custom_solver_factory; + void validate() const; }; @@ -143,11 +149,11 @@ class SolveBlock inline void SolveBlock::validate() const { - AssertThrow( - solve_type == SolveType::Constant || solve_type == SolveType::Explicit || - solve_type == SolveType::Linear || solve_type == SolveType::Newton, - dealii::ExcMessage( - "A valid solve type must be selected (Constant | Explicit | Linear | Newton)\n")); + AssertThrow(solve_type == SolveType::Constant || solve_type == SolveType::Explicit || + solve_type == SolveType::Linear || solve_type == SolveType::Newton || + solve_type == SolveType::CustomSolver, + dealii::ExcMessage("A valid solve type must be selected (Constant | " + "Explicit | Linear | Newton | CustomSolver)\n")); AssertThrow(!field_indices.empty(), dealii::ExcMessage("This solve block must manage at least 1 field.\n")); if (solve_type == SolveType::Newton) @@ -199,6 +205,14 @@ SolveBlock::validate() const dealii::ExcMessage("Constant \"solves\" do not have an RHS or LHS, " "and should have no dependencies.\n")); } + else if (solve_type == SolveType::CustomSolver) + { + AssertThrow(custom_solver_factory.has_value(), + dealii::ExcMessage( + "SolveType::CustomSolver requires setting custom_solver_factory.\n")); + AssertThrow(!field_indices.empty(), + dealii::ExcMessage("This solve block must manage at least 1 field.\n")); + } for (const auto &[field_index, dependency] : dependencies_rhs) { AssertThrow(dependency.src_flag == EvalFlags::nothing, diff --git a/include/prismspf/core/type_enums.h b/include/prismspf/core/type_enums.h index e9e59d24a..dee7cd4ad 100644 --- a/include/prismspf/core/type_enums.h +++ b/include/prismspf/core/type_enums.h @@ -40,7 +40,11 @@ enum SolveType : std::uint8_t * x -> x + damping * Deltax. * See [Newton's Method](https://en.wikipedia.org/wiki/Newton%27s_method) */ - Newton + Newton, + /** + * Custom solver to provide users greater flexibility + */ + CustomSolver }; /** diff --git a/include/prismspf/solvers/explicit_solver.h b/include/prismspf/solvers/explicit_solver.h index 0920f3959..6306f6793 100644 --- a/include/prismspf/solvers/explicit_solver.h +++ b/include/prismspf/solvers/explicit_solver.h @@ -44,9 +44,9 @@ class ExplicitSolver : public SolverBase {} void - init(const std::list &all_dependeny_sets) override + init(const std::list &all_dependency_sets) override { - SolverBase::init(all_dependeny_sets); + SolverBase::init(all_dependency_sets); unsigned int num_levels = solve_context->get_dof_manager().get_dof_handlers_levels().size(); // Initialize rhs_operators diff --git a/include/prismspf/solvers/linear_solver.h b/include/prismspf/solvers/linear_solver.h index 75fa0c00f..fdd71f372 100644 --- a/include/prismspf/solvers/linear_solver.h +++ b/include/prismspf/solvers/linear_solver.h @@ -71,9 +71,9 @@ class LinearSolver : public SolverBase * @brief Initialize the solver. */ void - init(const std::list &all_dependeny_sets) override + init(const std::list &all_dependency_sets) override { - SolverBase::init(all_dependeny_sets); + SolverBase::init(all_dependency_sets); rhs_vector.reinit(solutions.get_solution_full_vector(0)); // Initialize rhs_operator diff --git a/include/prismspf/solvers/newton_solver.h b/include/prismspf/solvers/newton_solver.h index bd99484a2..0d314e329 100644 --- a/include/prismspf/solvers/newton_solver.h +++ b/include/prismspf/solvers/newton_solver.h @@ -50,9 +50,9 @@ class NewtonSolver : public LinearSolver * @brief Initialize the solver. */ void - init(const std::list &all_dependeny_sets) override + init(const std::list &all_dependency_sets) override { - LinearSolver::init(all_dependeny_sets); + LinearSolver::init(all_dependency_sets); newton_update.reinit(solutions.get_solution_full_vector(0)); } diff --git a/include/prismspf/solvers/solver_base.h b/include/prismspf/solvers/solver_base.h index e07c203c4..ff3660d63 100644 --- a/include/prismspf/solvers/solver_base.h +++ b/include/prismspf/solvers/solver_base.h @@ -81,9 +81,9 @@ class SolverBase * @brief Initialize the solver. */ virtual void - init(const std::list &all_dependeny_sets) + init(const std::list &all_dependency_sets) { - DependencyExtents extents(solve_block.field_indices, all_dependeny_sets); + DependencyExtents extents(solve_block.field_indices, all_dependency_sets); solutions.init(extents.oldest_age); // Apply constraints. diff --git a/src/core/problem.cc b/src/core/problem.cc index 3ad81a272..e346168f3 100644 --- a/src/core/problem.cc +++ b/src/core/problem.cc @@ -57,6 +57,18 @@ make_solvers(const std::vector &solve_blocks, std::make_shared>(solve_block, solve_context)); break; + case SolveType::CustomSolver: + { + using FactoryFunc = + std::function>( + const SolveBlock &, + const SolveContext &)>; + + const auto &factory = + std::any_cast(solve_block.custom_solver_factory); + solvers.emplace_back(factory(solve_block, solve_context)); + break; + } default: AssertThrow(false, dealii::ExcMessage("Unknown solver type")); }