Skip to content

Commit bb7452a

Browse files
jijnjijn
committed
add application: allen_cahn_conserved
1 parent 325d973 commit bb7452a

5 files changed

Lines changed: 472 additions & 0 deletions

File tree

‎applications/allen_cahn_conserved/CMakeLists.txt‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
##
2+
# CMake script for the PRISMS-PF applications
3+
##
4+
5+
cmake_minimum_required(VERSION 3.28)
6+
cmake_policy(VERSION 3.28)
7+
8+
# Name the application
9+
set(APPLICATION_NAME "allen_cahn_conserved")
10+
11+
# Create a project for the application
12+
project(${APPLICATION_NAME} CXX)
13+
14+
# Set the build type to debug is not specified
15+
if(NOT CMAKE_BUILD_TYPE)
16+
set(CMAKE_BUILD_TYPE "Debug" CACHE STRING "Build type" FORCE)
17+
endif()
18+
19+
# Export compile commands for IDEs
20+
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
21+
22+
# Find the PRISMS-PF package
23+
find_package(
24+
PRISMS-PF
25+
REQUIRED
26+
HINTS
27+
${PRISMS_PF_DIR}
28+
$ENV{PRISMS_PF_DIR}
29+
)
30+
31+
# Declare the application executable
32+
add_executable(${APPLICATION_NAME})
33+
34+
# Rename the output executable to main
35+
set_target_properties(
36+
${APPLICATION_NAME}
37+
PROPERTIES
38+
OUTPUT_NAME
39+
main
40+
)
41+
42+
# Add sources
43+
target_sources(
44+
${APPLICATION_NAME}
45+
PRIVATE
46+
main.cc
47+
custom_pde.h
48+
)
49+
50+
# Require C++20
51+
# TODO: Remove this
52+
target_compile_features(${APPLICATION_NAME} PRIVATE cxx_std_20)
53+
set_target_properties(
54+
${APPLICATION_NAME}
55+
PROPERTIES
56+
CXX_STANDARD
57+
20
58+
CXX_STANDARD_REQUIRED
59+
ON
60+
CXX_EXTENSIONS
61+
OFF
62+
)
63+
64+
# Link PRISMS-PF
65+
target_link_libraries(
66+
${APPLICATION_NAME}
67+
PRIVATE
68+
$<$<CONFIG:Debug>:PRISMS-PF::prisms_pf_debug>
69+
$<$<CONFIG:Release>:PRISMS-PF::prisms_pf_release>
70+
)

‎applications/allen_cahn_conserved/allen_cahn_conserved.md‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
# PRISMS PhaseField: Globally Conserved Allen-Cahn Dynamics
2+
3+
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.
4+
5+
Consider a free energy expression of the form:
6+
7+
$$
8+
\begin{equation}
9+
F[c] = \int_{\Omega} f(c) + \frac{\kappa}{2} \nabla c \cdot \nabla c \,\mathrm{d} V,
10+
\end{equation}
11+
$$
12+
13+
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
14+
15+
$$
16+
\begin{equation}
17+
\int_{\Omega} c \,\mathrm{d} V = \mathrm{constant}.
18+
\end{equation}
19+
$$
20+
21+
The kinetics is given by the Allen-Cahn equation, which is the gradient flow of the functional $F[c]$:
22+
23+
$$
24+
\begin{equation}
25+
\frac{\partial c}{\partial t} = - M \mu,
26+
\end{equation}
27+
$$
28+
29+
where $M$ is the mobility and $\mu$ is the chemical potential defined as
30+
31+
$$
32+
\begin{equation}
33+
\mu = \frac{\delta F}{\delta c} = f_{,c}(c) - \kappa \nabla^2 c.
34+
\end{equation}
35+
$$
36+
37+
To conserve $c$, we can introduce a global Lagrange multiplier $\lambda(t)$
38+
39+
$$
40+
\begin{equation}
41+
\frac{\partial c}{\partial t} = - M (\mu - \lambda(t)).
42+
\end{equation}
43+
$$
44+
45+
The Lagrange multiplier can be solved from the conservation of $c$ (Eq. 2) as
46+
47+
$$
48+
\begin{equation}
49+
\lambda(t) = \frac{1}{V} \int_{\Omega} \mu \,\mathrm{d} V, \quad V = \int_{\Omega} \,\mathrm{d} V.
50+
\end{equation}
51+
$$
52+
53+
The governing equations are given by Eqs. 5, 4, and 6.
54+
55+
56+
## Time discretization
57+
Considering forward Euler explicit time stepping, we have the time discretized kinetics equations:
58+
59+
$$
60+
\begin{equation}
61+
c^{n+1} = c^{n} - \Delta t M(\mu^n-\lambda^n)
62+
\end{equation}
63+
$$
64+
65+
and
66+
67+
$$
68+
\begin{equation}
69+
\mu^{n} = f_{,c}^n - \kappa \nabla^2 c^n.
70+
\end{equation}
71+
$$
72+
73+
## Weak formulation
74+
In the weak formulation, considering an arbitrary variation $w$, the above equation can be expressed as a residual equations:
75+
76+
$$
77+
\begin{equation}
78+
\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
79+
\end{equation}
80+
$$
81+
82+
$$
83+
\begin{equation}
84+
r_{c} = c^{n} - \Delta t M(\mu^n-\lambda^n)
85+
\end{equation}
86+
$$
87+
88+
and
89+
90+
$$
91+
\begin{align}
92+
\int_{\Omega} w \mu^{n} \,\mathrm{d}V &= \int_{\Omega} w f_{,c}^n - w \kappa \nabla^2 c^{n} \,\mathrm{d}V \\
93+
&= \int_{\Omega} w ~ (f_{,c}^{n}) + \nabla w \cdot (\kappa \nabla c^{n}) \,\mathrm{d}V,
94+
\end{align}
95+
$$
96+
97+
$$
98+
\begin{equation}
99+
r_{\mu} = f_{,\eta}^{n}
100+
\end{equation}
101+
$$
102+
103+
$$
104+
\begin{equation}
105+
r_{\mu x} = \kappa \nabla \eta^{n}
106+
\end{equation}
107+
$$
108+
109+
where the Lagrange multiplier $\lambda^n$ for time step $n$ is calculated as
110+
111+
$$
112+
\begin{align}
113+
\lambda^n& = \frac{1}{V}\int_\Omega \mu^n \,\mathrm{d}V.
114+
\end{align}
115+
$$
116+
117+
The above values of $r_{c}$, $r_{\mu}$, and $r_{\mu x}$ are used in `compute_rhs()` in the following source file:
118+
`applications/allen_cahn_conserved/custom_pde.cc`

‎applications/allen_cahn_conserved/custom_pde.h‎

Lines changed: 167 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,167 @@
1+
// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
2+
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
3+
4+
#include <prismspf/core/pde_operator_base.h>
5+
6+
#include <prismspf/solvers/solver_base.h>
7+
8+
#include <prismspf/utilities/integrator.h>
9+
10+
#include <cmath>
11+
12+
PRISMS_PF_BEGIN_NAMESPACE
13+
14+
template <unsigned int dim, unsigned int degree, typename number>
15+
class CustomPDE : public PDEOperatorBase<dim, degree, number>
16+
{
17+
public:
18+
using ScalarValue = dealii::VectorizedArray<number>;
19+
using ScalarGrad = dealii::Tensor<1, dim, ScalarValue>;
20+
using ScalarHess = dealii::Tensor<2, dim, ScalarValue>;
21+
using VectorValue = dealii::Tensor<1, dim, ScalarValue>;
22+
using VectorGrad = dealii::Tensor<2, dim, ScalarValue>;
23+
using VectorHess = dealii::Tensor<3, dim, ScalarValue>;
24+
using PDEOperatorBase<dim, degree, number>::get_user_inputs;
25+
using PDEOperatorBase<dim, degree, number>::get_pf_tools;
26+
27+
/**
28+
* @brief Constructor.
29+
*/
30+
explicit CustomPDE(const UserInputParameters<dim> &_user_inputs,
31+
PhaseFieldTools<dim> &_pf_tools)
32+
: PDEOperatorBase<dim, degree, number>(_user_inputs, _pf_tools)
33+
, mobility(get_user_inputs().user_constants.get_double("mobility"))
34+
, kappa(get_user_inputs().user_constants.get_double("kappa"))
35+
{}
36+
37+
void
38+
set_lambda(number v)
39+
{
40+
lambda = v;
41+
}
42+
43+
private:
44+
number mobility;
45+
number kappa;
46+
number lambda;
47+
48+
void
49+
set_initial_condition([[maybe_unused]] const unsigned int &index,
50+
[[maybe_unused]] const unsigned int &component,
51+
[[maybe_unused]] const dealii::Point<dim> &point,
52+
[[maybe_unused]] number &scalar_value,
53+
[[maybe_unused]] number &vector_component_value) const override
54+
{
55+
const dealii::Tensor<1, dim> &mesh_size =
56+
get_user_inputs().spatial_discretization.rectangular_mesh.size;
57+
constexpr number center[12][3] = {
58+
{0.1, 0.3, 0},
59+
{0.8, 0.7, 0},
60+
{0.5, 0.2, 0},
61+
{0.4, 0.4, 0},
62+
{0.3, 0.9, 0},
63+
{0.8, 0.1, 0},
64+
{0.9, 0.5, 0},
65+
{0.0, 0.1, 0},
66+
{0.1, 0.6, 0},
67+
{0.5, 0.6, 0},
68+
{1, 1, 0},
69+
{0.7, 0.95, 0}
70+
};
71+
constexpr number rad[12] = {12, 14, 19, 16, 11, 12, 17, 15, 20, 10, 11, 14};
72+
73+
number dist = 0.0;
74+
number sdf = std::numeric_limits<double>::max();
75+
for (unsigned int i = 0; i < 12; i++)
76+
{
77+
dist = 0.0;
78+
for (unsigned int dir = 0; dir < dim; dir++)
79+
{
80+
const number comp_diff = point[dir] - center[i][dir] * mesh_size[dir];
81+
dist += comp_diff * comp_diff;
82+
}
83+
dist = std::sqrt(dist) - rad[i];
84+
sdf = std::min(sdf, dist);
85+
}
86+
scalar_value = 0.5 * (1.0 - std::tanh(sdf / 1.5));
87+
}
88+
89+
void
90+
compute_rhs(FieldContainer<dim, degree, number> &variable_list,
91+
const SimulationTimer &sim_timer,
92+
unsigned int solve_block_id) const override
93+
{
94+
if (solve_block_id == 0) // explicit c
95+
{
96+
ScalarValue c = variable_list.template get_value<Scalar, OldOne>(0);
97+
ScalarValue mu = variable_list.template get_value<Scalar, OldOne>(1);
98+
99+
ScalarValue eq_c = c - sim_timer.get_timestep() * mobility * (mu - lambda);
100+
101+
variable_list.set_value_term(0, eq_c);
102+
}
103+
else if (solve_block_id == 1) // mu
104+
{
105+
ScalarValue c = variable_list.template get_value<Scalar, Current>(0);
106+
ScalarGrad c_grad = variable_list.template get_gradient<Scalar, Current>(0);
107+
ScalarValue df_well = 4.0 * c * (c - 1.0) * (c - 0.5);
108+
109+
variable_list.set_value_term(1, df_well);
110+
variable_list.set_gradient_term(1, kappa * c_grad);
111+
}
112+
else if (solve_block_id == 2) // custom
113+
{
114+
}
115+
else if (solve_block_id == 3) // postprocess
116+
{
117+
ScalarValue c = variable_list.template get_value<Scalar, Current>(0);
118+
ScalarGrad c_grad = variable_list.template get_gradient<Scalar, Current>(0);
119+
ScalarValue f_chem = c * c * (1.0 - c) * (1.0 - c);
120+
ScalarValue f_grad = 0.5 * kappa * c_grad.norm_square();
121+
122+
variable_list.set_value_term(3, c_grad.norm());
123+
variable_list.set_value_term(4, f_chem + f_grad);
124+
}
125+
}
126+
};
127+
128+
template <unsigned int dim, unsigned int degree, typename number>
129+
class MyCustomSolver : public SolverBase<dim, degree, number>
130+
{
131+
private:
132+
CustomPDE<dim, degree, number> &custom_pde;
133+
134+
public:
135+
MyCustomSolver(const SolveBlock &solve_block,
136+
const SolveContext<dim, degree, number> &solve_context,
137+
CustomPDE<dim, degree, number> &_custom_pde)
138+
: SolverBase<dim, degree, number>(solve_block, solve_context)
139+
, custom_pde(_custom_pde)
140+
{}
141+
142+
void
143+
init(const std::list<DependencyMap> &all_dependency_sets) override
144+
{
145+
SolverBase<dim, degree, number>::init(all_dependency_sets);
146+
}
147+
148+
void
149+
solve_level(unsigned int relative_level) override
150+
{
151+
number mu_int = Integrator<dim, degree, number>::template integrate<0>(
152+
this->solve_context->get_dof_manager().get_field_dof_handler(1),
153+
this->solve_context->get_solution_indexer().get_solution_vector(1))[0];
154+
155+
const dealii::Tensor<1, dim> &mesh_size =
156+
this->solve_context->get_user_inputs().spatial_discretization.rectangular_mesh.size;
157+
number V = mesh_size[0] * mesh_size[1];
158+
159+
custom_pde.set_lambda(mu_int / V);
160+
}
161+
162+
void
163+
update() override
164+
{}
165+
};
166+
167+
PRISMS_PF_END_NAMESPACE

0 commit comments

Comments
 (0)