1 #ifndef INTEGRATOR_TOPOP_H
2 #define INTEGRATOR_TOPOP_H
9 #include "AMReX_ParallelDescriptor.H"
10 #include "AMReX_ParmParse.H"
27 #include "Numeric/Stencil.H"
33 #include "Operator/Operator.H"
61 if (type ==
"ellipse") value.
ic_psi =
new IC::Ellipse(value.geom, pp,
"psi.ic.ellipse");
62 else if (type ==
"constant") value.
ic_psi =
new IC::Constant(value.geom, pp,
"psi.ic.constant");
81 AMREX_D_TERM((value.geom[0].ProbHi()[0] - value.geom[0].ProbLo()[0]),
82 *(value.geom[0].ProbHi()[1] - value.geom[0].ProbLo()[1]),
83 *(value.geom[0].ProbHi()[2] - value.geom[0].ProbLo()[2]));
108 if (a_step > 0)
return;
110 for (
int lev = 0; lev <= finest_level; ++lev)
123 BL_PROFILE(
"TopOp::Advance");
127 amrex::Box domain = geom[lev].Domain();
132 for (amrex::MFIter mfi(*
psi_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
134 amrex::Box bx = mfi.tilebox();
137 amrex::Array4<const Set::Matrix>
const& sig = (*
stress_mf[lev]).array(mfi);
138 amrex::Array4<const Set::Matrix>
const& eps = (*
strain_mf[lev]).array(mfi);
139 amrex::Array4<const Set::Scalar>
const& psi = (*
psi_old_mf[lev]).array(mfi);
140 amrex::Array4<Set::Scalar>
const& psinew = (*
psi_mf[lev]).array(mfi);
142 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k)
146 driving_force +=
alpha * 2.0 * psi(i, j, k) * (2.0 * psi(i, j, k) * psi(i, j, k) - 3.0 * psi(i, j, k) + 1.0);
152 driving_force += -
gamma * 0.5 * (sig_avg.transpose() * eps_avg).trace() * psi(i, j, k);
156 psinew(i, j, k) = psi(i, j, k) - Lnow *
dt * driving_force;
157 if (psinew(i, j, k) < 0.0) psinew(i, j, k) = 0.0;
158 if (psinew(i, j, k) > 1.0) psinew(i, j, k) = 1.0;
164 const amrex::MFIter& mfi,
const amrex::Box&
box)
override
166 BL_PROFILE(
"TopOp::Integrate");
169 const amrex::Real* DX = geom[amrlev].CellSize();
170 Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
172 amrex::Array4<amrex::Real>
const& psi = (*
psi_mf[amrlev]).array(mfi);
173 amrex::Array4<const Set::Matrix>
const& sig = (*
stress_mf[amrlev]).array(mfi);
174 amrex::Array4<const Set::Matrix>
const& eps = (*
strain_mf[amrlev]).array(mfi);
175 amrex::ParallelFor(
box, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k)
177 volume += psi(i, j, k, 0) * dv;
178 w_chem_potential +=
alpha * psi(i, j, k, 0) * psi(i, j, k, 0) * (1. - psi(i, j, k, 0) * psi(i, j, k, 0)) * dv;
183 w_elastic +=
gamma * 0.5 * (sig_avg.transpose() * eps_avg).trace() * psi(i, j, k) * dv;
195 for (amrex::MFIter mfi(*
model_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
197 amrex::Box bx = mfi.tilebox();
198 amrex::Array4<char>
const& tags = a_tags.array(mfi);
201 amrex::Array4<Set::Scalar>
const& psi =
psi_mf[lev]->array(mfi);
202 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k)
208 tags(i, j, k) = amrex::TagBox::SET;