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