Line data Source code
1 : #ifndef INTEGRATOR_TOPOP_H
2 : #define INTEGRATOR_TOPOP_H
3 : #include <iostream>
4 : #include <fstream>
5 : #include <iomanip>
6 : #include <numeric>
7 :
8 : #include "AMReX.H"
9 : #include "AMReX_ParallelDescriptor.H"
10 : #include "AMReX_ParmParse.H"
11 :
12 : #include "IO/ParmParse.H"
13 : #include "Integrator/Base/Mechanics.H"
14 :
15 :
16 : #include "IC/IC.H"
17 : #include "BC/BC.H"
18 : #include "BC/Operator/Elastic/Constant.H"
19 : #include "BC/Operator/Elastic/TensionTest.H"
20 : #include "BC/Operator/Elastic/Expression.H"
21 :
22 : #include "IC/Ellipse.H"
23 : #include "IC/Voronoi.H"
24 : #include "IC/Constant.H"
25 : #include "IC/BMP.H"
26 : #include "BC/Constant.H"
27 : #include "Numeric/Stencil.H"
28 :
29 : #include "Model/Solid/Solid.H"
30 : #include "Solver/Nonlocal/Linear.H"
31 : #include "Solver/Nonlocal/Newton.H"
32 :
33 : #include "Operator/Operator.H"
34 :
35 :
36 : namespace Integrator
37 : {
38 : template<class MODEL>
39 : class TopOp: virtual public Base::Mechanics<MODEL>
40 : {
41 : public:
42 :
43 : TopOp(): Base::Mechanics<MODEL>() {}
44 1 : TopOp(IO::ParmParse& pp): Base::Mechanics<MODEL>()
45 : {
46 1 : Parse(*this, pp);
47 1 : }
48 :
49 : // The mechanics integrator manages the solution of an elastic
50 : // solve using the MLMG solver.
51 1 : static void Parse(TopOp& value, IO::ParmParse& pp)
52 : {
53 1 : Base::Mechanics<MODEL>::Parse(value, pp);
54 :
55 1 : pp_queryclass("model", value.model);
56 : // Read in IC for psi
57 1 : if (pp.contains("psi.ic.type"))
58 : {
59 1 : std::string type;
60 1 : pp_query_validate("psi.ic.type", type, {"ellipse", "constant"}); // Read IC type for the eta field
61 1 : if (type == "ellipse") value.ic_psi = new IC::Ellipse(value.geom, pp, "psi.ic.ellipse");
62 1 : else if (type == "constant") value.ic_psi = new IC::Constant(value.geom, pp, "psi.ic.constant");
63 0 : else Util::Abort(INFO, "Invalid value for psi.ic.type: ", type);
64 :
65 1 : value.bc_psi = new BC::Constant(1, pp, "psi.bc");
66 1 : value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, 2, "psi", true);
67 1 : value.RegisterNewFab(value.psi_old_mf, value.bc_psi, 1, 2, "psiold", false);
68 1 : value.psi_on = true;
69 : }
70 1 : pp_query_default("eta_ref_threshold", value.m_eta_ref_threshold, 0.01); // Refinement threshold based on eta
71 1 : pp_query_default("alpha", value.alpha, 1.0); // :math:`\alpha` parameter
72 1 : pp_query_default("beta", value.beta, 1.0); // :math:`\beta` parameter
73 1 : pp_query_default("gamma", value.gamma, 1.0); // :math:`\gamma` parameter
74 1 : pp_queryclass("L", value.L); // Mobility
75 1 : if (pp.contains("volume0frac"))
76 : {
77 0 : if (pp.contains("volume0")) Util::Abort(INFO, "Cannot specify volume0frac and volume0");
78 : Set::Scalar volumefrac;
79 0 : pp_query_default("volume0frac", volumefrac, 0.5); // Prescribed volume fraction
80 0 : value.volume0 = volumefrac *
81 0 : 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]));
84 : }
85 : else
86 1 : pp_query_default("volume0", value.volume0, 0.5); // Prescribed total vlume
87 :
88 1 : pp_queryclass("lambda", value.lambda); // Lagrange multiplier (can be interplated)
89 :
90 :
91 1 : value.RegisterIntegratedVariable(&value.volume, "volume");
92 1 : value.RegisterIntegratedVariable(&value.w_chem_potential, "chem_potential");
93 1 : value.RegisterIntegratedVariable(&value.w_bndry, "bndry");
94 1 : value.RegisterIntegratedVariable(&value.w_elastic, "elastic");
95 1 : }
96 :
97 1 : void Initialize(int lev) override
98 : {
99 1 : Base::Mechanics<MODEL>::Initialize(lev);
100 1 : if (psi_on) ic_psi->Initialize(lev, psi_mf);
101 1 : if (psi_on) ic_psi->Initialize(lev, psi_old_mf);
102 1 : }
103 :
104 2 : virtual void UpdateModel(int a_step, Set::Scalar /*a_time*/) override
105 : {
106 2 : if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
107 :
108 2 : if (a_step > 0) return;
109 :
110 2 : for (int lev = 0; lev <= finest_level; ++lev)
111 : {
112 1 : model_mf[lev]->setVal(model);
113 1 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
114 1 : Util::RealFillBoundary(*psi_mf[lev], geom[lev]);
115 1 : Util::RealFillBoundary(*psi_old_mf[lev], geom[lev]);
116 : }
117 :
118 : }
119 :
120 :
121 6 : void Advance(int lev, Set::Scalar time, Set::Scalar dt) override
122 : {
123 : BL_PROFILE("TopOp::Advance");
124 6 : Base::Mechanics<MODEL>::Advance(lev, time, dt);
125 6 : std::swap(psi_old_mf[lev], psi_mf[lev]);
126 6 : const Set::Scalar* DX = geom[lev].CellSize();
127 6 : amrex::Box domain = geom[lev].Domain();
128 :
129 6 : Set::Scalar Lnow = L(time);
130 6 : Set::Scalar lambdaT = lambda(time);
131 :
132 60 : for (amrex::MFIter mfi(*psi_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
133 : {
134 54 : amrex::Box bx = mfi.tilebox();
135 54 : bx.grow(1);
136 54 : bx = bx & domain;
137 54 : amrex::Array4<const Set::Matrix> const& sig = (*stress_mf[lev]).array(mfi);
138 54 : amrex::Array4<const Set::Matrix> const& eps = (*strain_mf[lev]).array(mfi);
139 54 : amrex::Array4<const Set::Scalar> const& psi = (*psi_old_mf[lev]).array(mfi);
140 54 : amrex::Array4<Set::Scalar> const& psinew = (*psi_mf[lev]).array(mfi);
141 :
142 12966 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
143 : {
144 12912 : Set::Scalar driving_force = 0.0;
145 :
146 51648 : 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);
147 12912 : driving_force += -beta * Numeric::Laplacian(psi, i, j, k, 0, DX);
148 :
149 12912 : Set::Matrix sig_avg = Numeric::Interpolate::NodeToCellAverage(sig, i, j, k, 0);
150 12912 : Set::Matrix eps_avg = Numeric::Interpolate::NodeToCellAverage(eps, i, j, k, 0);
151 :
152 12912 : driving_force += -gamma * 0.5 * (sig_avg.transpose() * eps_avg).trace() * psi(i, j, k);
153 :
154 12912 : driving_force += lambdaT * (volume - volume0);
155 :
156 25824 : psinew(i, j, k) = psi(i, j, k) - Lnow * dt * driving_force;
157 25824 : if (psinew(i, j, k) < 0.0) psinew(i, j, k) = 0.0;
158 28590 : if (psinew(i, j, k) > 1.0) psinew(i, j, k) = 1.0;
159 : });
160 : }
161 6 : }
162 :
163 2 : void Integrate(int amrlev, Set::Scalar time, int step,
164 : const amrex::MFIter& mfi, const amrex::Box& box) override
165 : {
166 : BL_PROFILE("TopOp::Integrate");
167 2 : Base::Mechanics<MODEL>::Integrate(amrlev, time, step, mfi, box);
168 :
169 2 : const amrex::Real* DX = geom[amrlev].CellSize();
170 2 : Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
171 :
172 2 : amrex::Array4<amrex::Real> const& psi = (*psi_mf[amrlev]).array(mfi);
173 2 : amrex::Array4<const Set::Matrix> const& sig = (*stress_mf[amrlev]).array(mfi);
174 2 : amrex::Array4<const Set::Matrix> const& eps = (*strain_mf[amrlev]).array(mfi);
175 2 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
176 : {
177 4096 : volume += psi(i, j, k, 0) * dv;
178 16384 : 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;
179 8192 : w_bndry += beta * 0.5 * Numeric::Gradient(psi, i, j, k, 0, DX).squaredNorm() * dv;
180 :
181 4096 : Set::Matrix sig_avg = Numeric::Interpolate::NodeToCellAverage(sig, i, j, k, 0);
182 4096 : Set::Matrix eps_avg = Numeric::Interpolate::NodeToCellAverage(eps, i, j, k, 0);
183 8192 : w_elastic += gamma * 0.5 * (sig_avg.transpose() * eps_avg).trace() * psi(i, j, k) * dv;
184 : });
185 2 : }
186 :
187 :
188 4 : void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar a_time, int a_ngrow) override
189 : {
190 4 : if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
191 4 : Base::Mechanics<MODEL>::TagCellsForRefinement(lev, a_tags, a_time, a_ngrow);
192 :
193 4 : Set::Vector DX(geom[lev].CellSize());
194 4 : Set::Scalar DXnorm = DX.lpNorm<2>();
195 17 : for (amrex::MFIter mfi(*model_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
196 : {
197 13 : amrex::Box bx = mfi.tilebox();
198 13 : amrex::Array4<char> const& tags = a_tags.array(mfi);
199 13 : if (psi_on)
200 : {
201 13 : amrex::Array4<Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
202 8218 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
203 : {
204 8205 : auto sten = Numeric::GetStencil(i, j, k, bx);
205 : {
206 8205 : Set::Vector gradpsi = Numeric::Gradient(psi, i, j, k, 0, DX.data(), sten);
207 8205 : if (gradpsi.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
208 1254 : tags(i, j, k) = amrex::TagBox::SET;
209 : }
210 : });
211 : }
212 : }
213 : }
214 :
215 :
216 : protected:
217 : MODEL model;
218 : IC::IC* ic_psi = nullptr;
219 : BC::BC<Set::Scalar>* bc_psi = nullptr;
220 : Set::Scalar m_eta_ref_threshold = NAN;
221 : Set::Field<Set::Scalar> psi_old_mf;
222 :
223 :
224 : using Base::Mechanics<MODEL>::m_type;
225 : using Base::Mechanics<MODEL>::finest_level;
226 : using Base::Mechanics<MODEL>::geom;
227 : using Base::Mechanics<MODEL>::model_mf;
228 : using Base::Mechanics<MODEL>::psi_mf;
229 : using Base::Mechanics<MODEL>::psi_on;
230 : using Base::Mechanics<MODEL>::stress_mf;
231 : using Base::Mechanics<MODEL>::strain_mf;
232 :
233 :
234 : Set::Scalar alpha = NAN;
235 : Set::Scalar beta = NAN;
236 : Set::Scalar gamma = NAN;
237 : //Set::Scalar L = 1.0;
238 : Numeric::Interpolator::Linear<Set::Scalar> L;
239 : Set::Scalar volume0 = NAN;
240 : //Set::Scalar lambda = 1.0;
241 : Numeric::Interpolator::Linear<Set::Scalar> lambda;
242 :
243 : Set::Scalar volume = NAN;
244 : Set::Scalar w_chem_potential = NAN;
245 : Set::Scalar w_bndry = NAN;
246 : Set::Scalar w_elastic = NAN;
247 : };
248 :
249 :
250 :
251 :
252 :
253 :
254 :
255 :
256 :
257 :
258 : } // namespace Integrator
259 : #endif
|
