Line data Source code
1 : #ifndef SOLVER_NONLOCAL_NEWTON
2 : #define SOLVER_NONLOCAL_NEWTON
3 :
4 : #include "Set/Set.H"
5 : #include "Operator/Elastic.H"
6 : #include "Solver/Nonlocal/Linear.H"
7 : #include "IO/ParmParse.H"
8 : #include "Numeric/Stencil.H"
9 :
10 : namespace Solver
11 : {
12 : namespace Nonlocal
13 : {
14 : template <typename T>
15 : class Newton: public Linear
16 : {
17 : public:
18 36 : Newton() {};
19 :
20 : Newton(Operator::Elastic<T::sym>& a_op)
21 : {
22 : this->Define(a_op);
23 : };
24 :
25 36 : ~Newton()
26 : {
27 36 : if (m_defined) Clear();
28 36 : }
29 :
30 734 : void Define(Operator::Elastic<T::sym>& a_op)
31 : {
32 734 : Linear::Define(a_op);
33 734 : m_elastic = dynamic_cast<Operator::Elastic<T::sym> *>(linop);
34 : //m_bc = &m_elastic->GetBC();
35 734 : }
36 734 : void Clear()
37 : {
38 734 : Linear::Clear();
39 734 : m_elastic = nullptr;
40 : //m_bc = nullptr;
41 734 : }
42 :
43 :
44 : void setNRIters(int a_nriters) { m_nriters = a_nriters; }
45 :
46 14 : void setPsi(Set::Field<Set::Scalar>& a_psi)
47 : {
48 14 : m_psi = &a_psi;
49 14 : }
50 :
51 : private:
52 : void prepareForSolve(const Set::Field<Set::Scalar>& a_u_mf,
53 : const Set::Field<Set::Scalar>& a_b_mf,
54 : Set::Field<Set::Scalar>& a_rhs_mf,
55 : Set::Field<Set::Matrix>& a_dw_mf,
56 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>>& a_ddw_mf,
57 : Set::Field<T>& a_model_mf)
58 : {
59 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
60 : {
61 : amrex::Box domain(linop->Geom(lev).growPeriodicDomain(1));
62 : domain.convert(amrex::IntVect::TheNodeVector());
63 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
64 : Set::Vector DX(linop->Geom(lev).CellSize());
65 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
66 : {
67 : amrex::Box bx = mfi.grownnodaltilebox();
68 : bx = bx & domain;
69 :
70 : amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
71 : amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
72 : amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
73 : amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, T::sym>> const& ddw = a_ddw_mf[lev]->array(mfi);
74 :
75 : // Set model internal dw and ddw.
76 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
77 : auto sten = Numeric::GetStencil(i, j, k, bx);
78 :
79 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
80 :
81 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
82 : {
83 : dw(i, j, k) = model(i, j, k).DW(gradu);
84 : ddw(i, j, k) = model(i, j, k).DDW(gradu);
85 : }
86 : else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
87 : {
88 : Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
89 : dw(i, j, k) = model(i, j, k).DW(eps);
90 : ddw(i, j, k) = model(i, j, k).DDW(eps);
91 : }
92 : else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
93 : {
94 : Set::Matrix F = gradu + Set::Matrix::Identity();
95 : dw(i, j, k) = model(i, j, k).DW(F);
96 : ddw(i, j, k) = model(i, j, k).DDW(F);
97 : }
98 : });
99 : }
100 :
101 : Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
102 : Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
103 : }
104 :
105 : m_elastic->SetModel(a_ddw_mf);
106 :
107 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
108 : {
109 : amrex::Box domain(linop->Geom(lev).growPeriodicDomain(1));
110 : domain.convert(amrex::IntVect::TheNodeVector());
111 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
112 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
113 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
114 : {
115 : amrex::Box bx = mfi.nodaltilebox();
116 : bx = bx & domain;
117 : amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
118 : amrex::Array4<const Set::Scalar> const& b = a_b_mf[lev]->array(mfi);
119 : amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
120 : amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
121 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
122 : {
123 : auto sten = Numeric::GetStencil(i, j, k, bx);
124 : // Do this if on the domain boundary
125 : if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
126 : {
127 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
128 : Set::Vector U(AMREX_D_DECL(u(i, j, k, 0), u(i, j, k, 1), u(i, j, k, 2)));
129 : Set::Vector ret = m_elastic->GetBC()(U, gradu, dw(i, j, k), i, j, k, bx);
130 : for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k, d) - ret(d);
131 : }
132 : else
133 : {
134 : Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx);
135 : for (int p = 0; p < AMREX_SPACEDIM; p++) rhs(i, j, k, p) = b(i, j, k, p) - divdw(p);
136 : }
137 : });
138 : }
139 : Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
140 : Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
141 : }
142 : }
143 :
144 1634 : void prepareForSolve(const Set::Field<Set::Vector>& a_u_mf,
145 : const Set::Field<Set::Vector>& a_b_mf,
146 : Set::Field<Set::Scalar>& a_rhs_mf,
147 : Set::Field<Set::Matrix>& a_dw_mf,
148 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>>& a_ddw_mf,
149 : Set::Field<T>& a_model_mf)
150 : {
151 3322 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
152 : {
153 1688 : amrex::Box domain(linop->Geom(lev).growPeriodicDomain(2));
154 1688 : domain.convert(amrex::IntVect::TheNodeVector());
155 1688 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
156 1688 : Set::Vector DX(linop->Geom(lev).CellSize());
157 3509 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
158 : {
159 1821 : amrex::Box bx = mfi.grownnodaltilebox();
160 1821 : bx = bx & domain;
161 :
162 1821 : amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
163 1821 : amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
164 1821 : amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
165 1821 : amrex::Array4<Set::Matrix4<AMREX_SPACEDIM, T::sym>> const& ddw = a_ddw_mf[lev]->array(mfi);
166 :
167 : // Set model internal dw and ddw.
168 545847 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
169 : {
170 276303 : auto sten = Numeric::GetStencil(i, j, k, bx);
171 :
172 276303 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
173 276303 : Set::Matrix kinvar;
174 552606 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
175 205269 : kinvar = gradu; // gradu
176 142068 : else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
177 0 : kinvar = 0.5 * (gradu + gradu.transpose()); // epsilon
178 142068 : else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
179 71034 : kinvar = gradu + Set::Matrix::Identity(); // F
180 :
181 828909 : dw(i, j, k) = model(i, j, k).DW(kinvar);
182 828909 : ddw(i, j, k) = model(i, j, k).DDW(kinvar);
183 :
184 : });
185 : }
186 :
187 1688 : a_dw_mf[lev]->setMultiGhost(true);
188 1688 : a_dw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
189 1688 : a_ddw_mf[lev]->setMultiGhost(true);
190 1688 : a_ddw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
191 : }
192 :
193 1634 : m_elastic->SetModel(a_ddw_mf);
194 1634 : if (m_psi)
195 58 : for (int i = 0; i <= a_b_mf.finest_level; i++)
196 44 : m_elastic->SetPsi(i, *(*m_psi)[i]);
197 :
198 :
199 3322 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
200 : {
201 1688 : amrex::Box domain(linop->Geom(lev).growPeriodicDomain(2));
202 1688 : domain.convert(amrex::IntVect::TheNodeVector());
203 1688 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
204 1688 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
205 3509 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
206 : {
207 1821 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
208 1821 : amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
209 1821 : amrex::Array4<const Set::Vector> const& b = a_b_mf[lev]->array(mfi);
210 1821 : amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
211 1821 : amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
212 :
213 : // This is for if psi is not being used and has not been set
214 1821 : if (m_psi == nullptr)
215 : {
216 189528 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
217 : {
218 187821 : auto sten = Numeric::GetStencil(i, j, k, bx);
219 : // Do this if on the domain boundary
220 187821 : if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
221 : {
222 77808 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
223 233424 : Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k), i, j, k, bx);
224 506640 : for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
225 77808 : }
226 : else
227 : {
228 110013 : Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
229 582465 : for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
230 : }
231 : });
232 : }
233 : // This is EXACTLY THE SAME AS THE ABOVE, excpet that we are now accounting
234 : // for psi terms. The reason we do this is because there's no good way (yet) to
235 : // initialize patches if psi = nullptr
236 : else
237 : {
238 114 : const amrex::Dim3 boxlo = amrex::lbound(bx), boxhi = amrex::ubound(bx);
239 114 : amrex::Array4<const Set::Scalar> const& psi = (*m_psi)[lev]->array(mfi);
240 88596 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
241 : {
242 88482 : auto sten = Numeric::GetStencil(i, j, k, bx);
243 : // Do this if on the domain boundary
244 88482 : if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
245 : {
246 2560 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
247 2560 : Set::Scalar psiavg = 1.0;
248 2560 : if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
249 : {
250 0 : psiavg = Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
251 : }
252 7680 : Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k) * psiavg, i, j, k, bx);
253 12800 : for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
254 2560 : }
255 : else
256 : {
257 85922 : Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
258 85922 : if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
259 : {
260 153860 : divdw *= Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
261 76930 : Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, dx);
262 153860 : divdw += dw(i, j, k) * gradpsi;
263 : }
264 429610 : for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
265 : }
266 : });
267 : }
268 : }
269 :
270 1688 : a_rhs_mf[lev]->setMultiGhost(true);
271 1688 : a_rhs_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
272 1688 : a_ddw_mf[lev]->setMultiGhost(true);
273 1688 : a_ddw_mf[lev]->FillBoundaryAndSync(m_elastic->Geom(lev).periodicity());
274 : }
275 1634 : }
276 :
277 :
278 : public:
279 734 : Set::Scalar solve(const Set::Field<Set::Vector>& a_u_mf,
280 : const Set::Field<Set::Vector>& a_b_mf,
281 : Set::Field<T>& a_model_mf,
282 : Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
283 : {
284 734 : Set::Field<Set::Scalar> dsol_mf, rhs_mf;
285 734 : Set::Field<Set::Matrix> dw_mf;
286 734 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
287 :
288 734 : dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
289 734 : dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
290 734 : ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
291 734 : rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
292 1522 : for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
293 : {
294 1576 : dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
295 788 : a_u_mf[lev]->DistributionMap(),
296 : a_u_mf.NComp(),
297 788 : a_u_mf[lev]->nGrow());
298 1576 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
299 788 : a_b_mf[lev]->DistributionMap(),
300 : 1,
301 788 : a_b_mf[lev]->nGrow());
302 1576 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
303 788 : a_b_mf[lev]->DistributionMap(),
304 : 1,
305 788 : a_b_mf[lev]->nGrow());
306 1576 : rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
307 788 : a_b_mf[lev]->DistributionMap(),
308 : a_b_mf.NComp(),
309 788 : a_b_mf[lev]->nGrow());
310 :
311 788 : dsol_mf[lev]->setVal(0.0);
312 788 : dw_mf[lev]->setVal(Set::Matrix::Zero());
313 788 : ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
314 :
315 788 : a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
316 : //amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
317 : }
318 :
319 2366 : for (int nriter = 0; nriter < m_nriters; nriter++)
320 : {
321 8170 : if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
322 :
323 1634 : prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
324 :
325 :
326 1634 : if (nriter == m_nriters) break;
327 :
328 1634 : Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
329 :
330 1634 : Set::Scalar cornorm = 0, solnorm = 0;
331 3322 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
332 : {
333 5664 : for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
334 : {
335 3976 : Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
336 3976 : if (tmpcornorm > cornorm) cornorm = tmpcornorm;
337 :
338 : //Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp,0);
339 : //if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
340 : }
341 :
342 : }
343 : Set::Scalar relnorm;
344 1634 : if (solnorm == 0) relnorm = cornorm;
345 0 : else relnorm = cornorm / solnorm;
346 8170 : if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
347 :
348 3322 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
349 1688 : a_u_mf.AddFrom(lev, *dsol_mf[lev], 0, 2);
350 : //amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
351 :
352 1634 : if (relnorm < m_nrtolerance)
353 2 : return relnorm;
354 :
355 : }
356 :
357 732 : return 0.0;
358 734 : }
359 :
360 : Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
361 : const Set::Field<Set::Scalar>& a_b_mf,
362 : Set::Field<T>& a_model_mf,
363 : Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
364 : {
365 : Set::Field<Set::Scalar> dsol_mf, rhs_mf;
366 : Set::Field<Set::Matrix> dw_mf;
367 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
368 :
369 : dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
370 : dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
371 : ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
372 : rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
373 : for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
374 : {
375 : dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
376 : a_u_mf[lev]->DistributionMap(),
377 : a_u_mf[lev]->nComp(),
378 : a_u_mf[lev]->nGrow());
379 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
380 : a_b_mf[lev]->DistributionMap(),
381 : 1,
382 : a_b_mf[lev]->nGrow());
383 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
384 : a_b_mf[lev]->DistributionMap(),
385 : 1,
386 : a_b_mf[lev]->nGrow());
387 : rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
388 : a_b_mf[lev]->DistributionMap(),
389 : a_b_mf[lev]->nComp(),
390 : a_b_mf[lev]->nGrow());
391 :
392 : dsol_mf[lev]->setVal(0.0);
393 : dw_mf[lev]->setVal(Set::Matrix::Zero());
394 : ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
395 :
396 : amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
397 : }
398 :
399 : for (int nriter = 0; nriter < m_nriters; nriter++)
400 : {
401 : if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
402 :
403 : prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
404 :
405 :
406 : if (nriter == m_nriters) break;
407 :
408 : Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
409 :
410 : Set::Scalar cornorm = 0, solnorm = 0;
411 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
412 : {
413 : for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
414 : {
415 : Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
416 : if (tmpcornorm > cornorm) cornorm = tmpcornorm;
417 :
418 : Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp, 0);
419 : if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
420 : }
421 :
422 : }
423 : Set::Scalar relnorm;
424 : if (solnorm == 0) relnorm = cornorm;
425 : else relnorm = cornorm / solnorm;
426 : if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
427 :
428 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
429 : amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
430 :
431 : if (relnorm < m_nrtolerance)
432 : return relnorm;
433 :
434 : }
435 :
436 : return 0.0;
437 : }
438 : Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
439 : const Set::Field<Set::Scalar>& a_b_mf,
440 : Set::Field<T>& a_model_mf)
441 : {
442 : return solve(a_u_mf, a_b_mf, a_model_mf, tol_rel, tol_abs);
443 : }
444 :
445 : void compResidual(Set::Field<Set::Scalar>& a_res_mf,
446 : Set::Field<Set::Scalar>& a_u_mf,
447 : Set::Field<Set::Scalar>& a_b_mf,
448 : Set::Field<T>& a_model_mf)
449 : {
450 : Set::Field<Set::Matrix> dw_mf;
451 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
452 : dw_mf.resize(a_u_mf.size());
453 : ddw_mf.resize(a_u_mf.size());
454 : for (int lev = 0; lev < a_u_mf.size(); lev++)
455 : {
456 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
457 : a_b_mf[lev]->DistributionMap(),
458 : 1, a_b_mf[lev]->nGrow());
459 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
460 : a_b_mf[lev]->DistributionMap(),
461 : 1, a_b_mf[lev]->nGrow());
462 : dw_mf[lev]->setVal(Set::Matrix::Zero());
463 : }
464 :
465 : prepareForSolve(a_u_mf, a_b_mf, a_res_mf, dw_mf, ddw_mf, a_model_mf);
466 : }
467 :
468 : void compResidual(Set::Field<Set::Vector>& a_res_mf,
469 : Set::Field<Set::Vector>& a_u_mf,
470 : Set::Field<Set::Vector>& a_b_mf,
471 : Set::Field<T>& a_model_mf)
472 : {
473 : Set::Field<Set::Matrix> dw_mf;
474 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
475 : Set::Field<Set::Scalar> res_mf(a_res_mf.size());
476 : dw_mf.resize(a_u_mf.size());
477 : ddw_mf.resize(a_u_mf.size());
478 : res_mf.resize(a_u_mf.size());
479 : for (int lev = 0; lev < a_u_mf.size(); lev++)
480 : {
481 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
482 : a_b_mf[lev]->DistributionMap(),
483 : 1, a_b_mf[lev]->nGrow());
484 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
485 : a_b_mf[lev]->DistributionMap(),
486 : 1, a_b_mf[lev]->nGrow());
487 : res_mf.Define(lev, a_b_mf[lev]->boxArray(),
488 : a_b_mf[lev]->DistributionMap(),
489 : AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
490 : dw_mf[lev]->setVal(Set::Matrix::Zero());
491 : }
492 :
493 : prepareForSolve(a_u_mf, a_b_mf, res_mf, dw_mf, ddw_mf, a_model_mf);
494 :
495 : for (int lev = 0; lev < a_res_mf.size(); ++lev)
496 : {
497 : Util::RealFillBoundary(*res_mf[lev], m_elastic->Geom(lev));
498 : a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
499 : }
500 : }
501 :
502 12 : void compLinearSolverResidual(Set::Field<Set::Vector>& a_res_mf,
503 : Set::Field<Set::Vector>& a_u_mf,
504 : Set::Field<Set::Vector>& a_b_mf)
505 : {
506 84 : Util::Assert(INFO, TEST(a_res_mf.finest_level == a_u_mf.finest_level));
507 84 : Util::Assert(INFO, TEST(a_u_mf.finest_level == a_b_mf.finest_level));
508 12 : Set::Field<Set::Scalar> res_mf(a_res_mf.finest_level + 1), sol_mf(a_u_mf.finest_level + 1), rhs_mf(a_b_mf.finest_level + 1);
509 42 : for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
510 : {
511 30 : res_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
512 30 : sol_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
513 30 : rhs_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
514 :
515 30 : a_u_mf.Copy(lev, *sol_mf[lev], 0, 2);
516 30 : a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
517 : }
518 :
519 12 : mlmg->compResidual(amrex::GetVecOfPtrs(res_mf), amrex::GetVecOfPtrs(sol_mf), amrex::GetVecOfConstPtrs(rhs_mf));
520 :
521 42 : for (int lev = 0; lev <= a_res_mf.finest_level; ++lev)
522 : {
523 30 : a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
524 : }
525 12 : }
526 :
527 :
528 : void W(Set::Field<Set::Scalar>& a_w_mf,
529 : Set::Field<Set::Scalar>& a_u_mf,
530 : Set::Field<T>& a_model_mf)
531 : {
532 : for (int lev = 0; lev < a_u_mf.size(); lev++)
533 : {
534 : BL_PROFILE("Solver::Nonlocal::Newton::DW()");
535 :
536 : const amrex::Real* DX = linop->Geom(lev).CellSize();
537 : amrex::Box domain(linop->Geom(lev).Domain());
538 : domain.convert(amrex::IntVect::TheNodeVector());
539 :
540 : for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
541 : {
542 : const Box& bx = mfi.tilebox();
543 : amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
544 : amrex::Array4<amrex::Real> const& w = a_w_mf[lev]->array(mfi);
545 : amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
546 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
547 : {
548 : Set::Matrix gradu;
549 :
550 : auto sten = Numeric::GetStencil(i, j, k, bx);
551 :
552 : // Fill gradu
553 : for (int p = 0; p < AMREX_SPACEDIM; p++)
554 : {
555 : AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
556 : gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
557 : gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
558 : }
559 :
560 : if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
561 : w(i, j, k) = C(i, j, k).W(gradu);
562 : else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
563 : w(i, j, k) = C(i, j, k).W(0.5 * (gradu + gradu.transpose()));
564 : else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
565 : w(i, j, k) = C(i, j, k).W(gradu + Set::Matrix::Identity());
566 : });
567 : }
568 : }
569 : }
570 :
571 : void DW(Set::Field<Set::Scalar>& a_dw_mf,
572 : Set::Field<Set::Scalar>& a_u_mf,
573 : Set::Field<T>& a_model_mf)
574 : {
575 : for (int lev = 0; lev < a_u_mf.size(); lev++)
576 : {
577 : BL_PROFILE("Solver::Nonlocal::Newton::DW()");
578 :
579 : const amrex::Real* DX = linop->Geom(lev).CellSize();
580 : amrex::Box domain(linop->Geom(lev).Domain());
581 : domain.convert(amrex::IntVect::TheNodeVector());
582 :
583 : for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
584 : {
585 : const Box& bx = mfi.tilebox();
586 : amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
587 : amrex::Array4<amrex::Real> const& dw = a_dw_mf[lev]->array(mfi);
588 : amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
589 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
590 : {
591 : Set::Matrix gradu;
592 :
593 : auto sten = Numeric::GetStencil(i, j, k, bx);
594 :
595 : // Fill gradu
596 : for (int p = 0; p < AMREX_SPACEDIM; p++)
597 : {
598 : AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
599 : gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
600 : gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
601 : }
602 :
603 : Set::Matrix sig = Set::Matrix::Zero();
604 :
605 : if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
606 : sig = C(i, j, k).DW(gradu);
607 : else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
608 : sig = C(i, j, k).DW(0.5 * (gradu + gradu.transpose()));
609 : else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
610 : sig = C(i, j, k).DW(gradu + Set::Matrix::Identity());
611 :
612 : // = C(i,j,k)(gradu,m_homogeneous);
613 :
614 : AMREX_D_PICK(dw(i, j, k, 0) = sig(0, 0);
615 : ,
616 : dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1);
617 : dw(i, j, k, 2) = sig(1, 0); dw(i, j, k, 3) = sig(1, 1);
618 : ,
619 : dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1); dw(i, j, k, 2) = sig(0, 2);
620 : dw(i, j, k, 3) = sig(1, 0); dw(i, j, k, 4) = sig(1, 1); dw(i, j, k, 5) = sig(1, 2);
621 : dw(i, j, k, 6) = sig(2, 0); dw(i, j, k, 7) = sig(2, 1); dw(i, j, k, 8) = sig(2, 2););
622 :
623 : });
624 : }
625 : }
626 : }
627 :
628 :
629 : public:
630 : int m_nriters = 1;
631 : Set::Scalar m_nrtolerance = 0.0;
632 : Operator::Elastic<T::sym>* m_elastic;
633 : //BC::Operator::Elastic::Elastic *m_bc;
634 :
635 : Set::Field<Set::Scalar>* m_psi = nullptr;
636 :
637 : public:
638 : // These paramters control a standard Newton-Raphson solve.
639 : //
640 : // **Note**:
641 : // This class inherits all of the linear solve paramters
642 : // from its parent class, :ref:`Solver::Nonlocal::Linear`
643 30 : static void Parse(Newton<T>& value, amrex::ParmParse& pp)
644 : {
645 30 : Linear::Parse(value, pp);
646 :
647 : // Number of newton-raphson iterations.
648 30 : pp_query("nriters", value.m_nriters);
649 :
650 : // Tolerance to use for newton-raphson convergence
651 30 : pp_query("nrtolerance", value.m_nrtolerance);
652 30 : }
653 :
654 : };
655 : } // namespace Nonlocal
656 : } // namespace Solver
657 :
658 :
659 : #endif
|
