Alamo
Newton.H
Go to the documentation of this file.
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
10namespace Solver
11{
12namespace Nonlocal
13{
14template <typename T>
15class Newton: public Linear
16{
17public:
18 Newton() {};
19
20 Newton(Operator::Elastic<T::sym>& a_op)
21 {
22 this->Define(a_op);
23 };
24
25 ~Newton()
26 {
27 if (m_defined) Clear();
28 }
29
30 void Define(Operator::Elastic<T::sym>& a_op)
31 {
32 Linear::Define(a_op);
33 m_elastic = dynamic_cast<Operator::Elastic<T::sym> *>(linop);
34 //m_bc = &m_elastic->GetBC();
35 }
36 void Clear()
37 {
38 Linear::Clear();
39 m_elastic = nullptr;
40 //m_bc = nullptr;
41 }
42
43
44 void setNRIters(int a_nriters) { m_nriters = a_nriters; }
45
46 void setPsi(Set::Field<Set::Scalar>& a_psi)
47 {
48 m_psi = &a_psi;
49 }
50
51private:
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).Domain());
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).Domain());
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 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 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
152 {
153 amrex::Box domain(linop->Geom(lev).Domain());
154 domain.convert(amrex::IntVect::TheNodeVector());
155 const Set::Scalar* dx = linop->Geom(lev).CellSize();
156 Set::Vector DX(linop->Geom(lev).CellSize());
157 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
158 {
159 amrex::Box bx = mfi.grownnodaltilebox();
160 bx = bx & domain;
161
162 amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
163 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
164 amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
165 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 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
169 {
170 auto sten = Numeric::GetStencil(i, j, k, bx);
171
172 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
173 Set::Matrix kinvar;
174 if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
175 kinvar = gradu; // gradu
176 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
177 kinvar = 0.5 * (gradu + gradu.transpose()); // epsilon
178 else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
179 kinvar = gradu + Set::Matrix::Identity(); // F
180
181 dw(i, j, k) = model(i, j, k).DW(kinvar);
182 ddw(i, j, k) = model(i, j, k).DDW(kinvar);
183
184 });
185 }
186
187 Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
188 Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
189 }
190
191 m_elastic->SetModel(a_ddw_mf);
192 if (m_psi)
193 for (int i = 0; i <= a_b_mf.finest_level; i++)
194 m_elastic->SetPsi(i, *(*m_psi)[i]);
195
196
197 for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
198 {
199 amrex::Box domain(linop->Geom(lev).Domain());
200 domain.convert(amrex::IntVect::TheNodeVector());
201 const Set::Scalar* dx = linop->Geom(lev).CellSize();
202 const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
203 for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
204 {
205 amrex::Box bx = mfi.grownnodaltilebox() & domain;
206 amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
207 amrex::Array4<const Set::Vector> const& b = a_b_mf[lev]->array(mfi);
208 amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
209 amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
210
211 // This is for if psi is not being used and has not been set
212 if (m_psi == nullptr)
213 {
214 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
215 {
216 auto sten = Numeric::GetStencil(i, j, k, bx);
217 // Do this if on the domain boundary
218 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
219 {
220 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
221 Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k), i, j, k, bx);
222 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
223 }
224 else
225 {
226 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
227 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
228 }
229 });
230 }
231 // This is EXACTLY THE SAME AS THE ABOVE, excpet that we are now accounting
232 // for psi terms. The reason we do this is because there's no good way (yet) to
233 // initialize patches if psi = nullptr
234 else
235 {
236 const amrex::Dim3 boxlo = amrex::lbound(bx), boxhi = amrex::ubound(bx);
237 amrex::Array4<const Set::Scalar> const& psi = (*m_psi)[lev]->array(mfi);
238 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
239 {
240 auto sten = Numeric::GetStencil(i, j, k, bx);
241 // Do this if on the domain boundary
242 if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
243 {
244 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
245 Set::Scalar psiavg = 1.0;
246 if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
247 {
248 psiavg = Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
249 }
250 Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k) * psiavg, i, j, k, bx);
251 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - ret(d);
252 }
253 else
254 {
255 Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
256 if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
257 {
258 divdw *= Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
259 Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, dx);
260 divdw += dw(i, j, k) * gradpsi;
261 }
262 for (int d = 0; d < AMREX_SPACEDIM; d++) rhs(i, j, k, d) = b(i, j, k)(d) - divdw(d);
263 }
264 });
265 }
266 }
267
268 Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
269 Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
270 }
271 }
272
273
274public:
275 Set::Scalar solve(const Set::Field<Set::Vector>& a_u_mf,
276 const Set::Field<Set::Vector>& a_b_mf,
277 Set::Field<T>& a_model_mf,
278 Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
279 {
280 Set::Field<Set::Scalar> dsol_mf, rhs_mf;
281 Set::Field<Set::Matrix> dw_mf;
282 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
283
284 dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
285 dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
286 ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
287 rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
288 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
289 {
290 dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
291 a_u_mf[lev]->DistributionMap(),
292 a_u_mf.NComp(),
293 a_u_mf[lev]->nGrow());
294 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
295 a_b_mf[lev]->DistributionMap(),
296 1,
297 a_b_mf[lev]->nGrow());
298 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
299 a_b_mf[lev]->DistributionMap(),
300 1,
301 a_b_mf[lev]->nGrow());
302 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
303 a_b_mf[lev]->DistributionMap(),
304 a_b_mf.NComp(),
305 a_b_mf[lev]->nGrow());
306
307 dsol_mf[lev]->setVal(0.0);
308 dw_mf[lev]->setVal(Set::Matrix::Zero());
309 ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
310
311 a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
312 //amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
313 }
314
315 for (int nriter = 0; nriter < m_nriters; nriter++)
316 {
317 if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
318
319 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
320
321
322 if (nriter == m_nriters) break;
323
324 Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
325
326 Set::Scalar cornorm = 0, solnorm = 0;
327 for (int lev = 0; lev < dsol_mf.size(); ++lev)
328 {
329 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
330 {
331 Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
332 if (tmpcornorm > cornorm) cornorm = tmpcornorm;
333
334 //Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp,0);
335 //if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
336 }
337
338 }
339 Set::Scalar relnorm;
340 if (solnorm == 0) relnorm = cornorm;
341 else relnorm = cornorm / solnorm;
342 if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
343
344 for (int lev = 0; lev < dsol_mf.size(); ++lev)
345 a_u_mf.AddFrom(lev, *dsol_mf[lev], 0, 2);
346 //amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
347
348 if (relnorm < m_nrtolerance)
349 return relnorm;
350
351 }
352
353 return 0.0;
354 }
355
356 Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
357 const Set::Field<Set::Scalar>& a_b_mf,
358 Set::Field<T>& a_model_mf,
359 Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr)
360 {
361 Set::Field<Set::Scalar> dsol_mf, rhs_mf;
362 Set::Field<Set::Matrix> dw_mf;
363 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
364
365 dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
366 dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
367 ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
368 rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
369 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
370 {
371 dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
372 a_u_mf[lev]->DistributionMap(),
373 a_u_mf[lev]->nComp(),
374 a_u_mf[lev]->nGrow());
375 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
376 a_b_mf[lev]->DistributionMap(),
377 1,
378 a_b_mf[lev]->nGrow());
379 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
380 a_b_mf[lev]->DistributionMap(),
381 1,
382 a_b_mf[lev]->nGrow());
383 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
384 a_b_mf[lev]->DistributionMap(),
385 a_b_mf[lev]->nComp(),
386 a_b_mf[lev]->nGrow());
387
388 dsol_mf[lev]->setVal(0.0);
389 dw_mf[lev]->setVal(Set::Matrix::Zero());
390 ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
391
392 amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
393 }
394
395 for (int nriter = 0; nriter < m_nriters; nriter++)
396 {
397 if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
398
399 prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
400
401
402 if (nriter == m_nriters) break;
403
404 Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
405 //Solver::Nonlocal::Linear::solve(GetVecOfPtrs(dsol_mf), GetVecOfConstPtrs(rhs_mf), a_tol_rel, a_tol_abs,checkpoint_file);
406
407 Set::Scalar cornorm = 0, solnorm = 0;
408 for (int lev = 0; lev < dsol_mf.size(); ++lev)
409 {
410 for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
411 {
412 Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
413 if (tmpcornorm > cornorm) cornorm = tmpcornorm;
414
415 Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp, 0);
416 if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
417 }
418
419 }
420 Set::Scalar relnorm;
421 if (solnorm == 0) relnorm = cornorm;
422 else relnorm = cornorm / solnorm;
423 if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
424
425 for (int lev = 0; lev < dsol_mf.size(); ++lev)
426 amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
427
428 if (relnorm < m_nrtolerance)
429 return relnorm;
430
431 }
432
433 return 0.0;
434 }
435 Set::Scalar solve(const Set::Field<Set::Scalar>& a_u_mf,
436 const Set::Field<Set::Scalar>& a_b_mf,
437 Set::Field<T>& a_model_mf)
438 {
439 return solve(a_u_mf, a_b_mf, a_model_mf, tol_rel, tol_abs);
440 }
441
442 void compResidual(Set::Field<Set::Scalar>& a_res_mf,
443 Set::Field<Set::Scalar>& a_u_mf,
444 Set::Field<Set::Scalar>& a_b_mf,
445 Set::Field<T>& a_model_mf)
446 {
447 Set::Field<Set::Matrix> dw_mf;
448 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
449 dw_mf.resize(a_u_mf.size());
450 ddw_mf.resize(a_u_mf.size());
451 for (int lev = 0; lev < a_u_mf.size(); lev++)
452 {
453 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
454 a_b_mf[lev]->DistributionMap(),
455 1, a_b_mf[lev]->nGrow());
456 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
457 a_b_mf[lev]->DistributionMap(),
458 1, a_b_mf[lev]->nGrow());
459 dw_mf[lev]->setVal(Set::Matrix::Zero());
460 }
461
462 prepareForSolve(a_u_mf, a_b_mf, a_res_mf, dw_mf, ddw_mf, a_model_mf);
463 }
464
465 void compResidual(Set::Field<Set::Vector>& a_res_mf,
466 Set::Field<Set::Vector>& a_u_mf,
467 Set::Field<Set::Vector>& a_b_mf,
468 Set::Field<T>& a_model_mf)
469 {
470 Set::Field<Set::Matrix> dw_mf;
471 Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
472 Set::Field<Set::Scalar> res_mf(a_res_mf.size());
473 dw_mf.resize(a_u_mf.size());
474 ddw_mf.resize(a_u_mf.size());
475 res_mf.resize(a_u_mf.size());
476 for (int lev = 0; lev < a_u_mf.size(); lev++)
477 {
478 dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
479 a_b_mf[lev]->DistributionMap(),
480 1, a_b_mf[lev]->nGrow());
481 ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
482 a_b_mf[lev]->DistributionMap(),
483 1, a_b_mf[lev]->nGrow());
484 res_mf.Define(lev, a_b_mf[lev]->boxArray(),
485 a_b_mf[lev]->DistributionMap(),
486 AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
487 dw_mf[lev]->setVal(Set::Matrix::Zero());
488 }
489
490 prepareForSolve(a_u_mf, a_b_mf, res_mf, dw_mf, ddw_mf, a_model_mf);
491
492 for (int lev = 0; lev < a_res_mf.size(); ++lev)
493 {
494 Util::RealFillBoundary(*res_mf[lev], m_elastic->Geom(lev));
495 a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
496 }
497 }
498
499 void compLinearSolverResidual(Set::Field<Set::Vector>& a_res_mf,
500 Set::Field<Set::Vector>& a_u_mf,
501 Set::Field<Set::Vector>& a_b_mf)
502 {
503 Util::Assert(INFO, TEST(a_res_mf.finest_level == a_u_mf.finest_level));
504 Util::Assert(INFO, TEST(a_u_mf.finest_level == a_b_mf.finest_level));
505 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);
506 for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
507 {
508 res_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
509 sol_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
510 rhs_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
511
512 a_u_mf.Copy(lev, *sol_mf[lev], 0, 2);
513 a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
514 }
515
516 mlmg->compResidual(amrex::GetVecOfPtrs(res_mf), amrex::GetVecOfPtrs(sol_mf), amrex::GetVecOfConstPtrs(rhs_mf));
517
518 for (int lev = 0; lev <= a_res_mf.finest_level; ++lev)
519 {
520 a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
521 }
522 }
523
524
525 void W(Set::Field<Set::Scalar>& a_w_mf,
526 Set::Field<Set::Scalar>& a_u_mf,
527 Set::Field<T>& a_model_mf)
528 {
529 for (int lev = 0; lev < a_u_mf.size(); lev++)
530 {
531 BL_PROFILE("Solver::Nonlocal::Newton::DW()");
532
533 const amrex::Real* DX = linop->Geom(lev).CellSize();
534 amrex::Box domain(linop->Geom(lev).Domain());
535 domain.convert(amrex::IntVect::TheNodeVector());
536
537 for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
538 {
539 const Box& bx = mfi.tilebox();
540 amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
541 amrex::Array4<amrex::Real> const& w = a_w_mf[lev]->array(mfi);
542 amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
543 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
544 {
545 Set::Matrix gradu;
546
547 auto sten = Numeric::GetStencil(i, j, k, bx);
548
549 // Fill gradu
550 for (int p = 0; p < AMREX_SPACEDIM; p++)
551 {
552 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
553 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
554 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
555 }
556
557 if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
558 w(i, j, k) = C(i, j, k).W(gradu);
559 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
560 w(i, j, k) = C(i, j, k).W(0.5 * (gradu + gradu.transpose()));
561 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
562 w(i, j, k) = C(i, j, k).W(gradu + Set::Matrix::Identity());
563 });
564 }
565 }
566 }
567
568 void DW(Set::Field<Set::Scalar>& a_dw_mf,
569 Set::Field<Set::Scalar>& a_u_mf,
570 Set::Field<T>& a_model_mf)
571 {
572 for (int lev = 0; lev < a_u_mf.size(); lev++)
573 {
574 BL_PROFILE("Solver::Nonlocal::Newton::DW()");
575
576 const amrex::Real* DX = linop->Geom(lev).CellSize();
577 amrex::Box domain(linop->Geom(lev).Domain());
578 domain.convert(amrex::IntVect::TheNodeVector());
579
580 for (MFIter mfi(*a_u_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
581 {
582 const Box& bx = mfi.tilebox();
583 amrex::Array4<T> const& C = a_model_mf[lev]->array(mfi);
584 amrex::Array4<amrex::Real> const& dw = a_dw_mf[lev]->array(mfi);
585 amrex::Array4<const amrex::Real> const& u = a_u_mf[lev]->array(mfi);
586 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
587 {
588 Set::Matrix gradu;
589
590 auto sten = Numeric::GetStencil(i, j, k, bx);
591
592 // Fill gradu
593 for (int p = 0; p < AMREX_SPACEDIM; p++)
594 {
595 AMREX_D_TERM(gradu(p, 0) = (Numeric::Stencil<Set::Scalar, 1, 0, 0>::D(u, i, j, k, p, DX, sten));,
596 gradu(p, 1) = (Numeric::Stencil<Set::Scalar, 0, 1, 0>::D(u, i, j, k, p, DX, sten));,
597 gradu(p, 2) = (Numeric::Stencil<Set::Scalar, 0, 0, 1>::D(u, i, j, k, p, DX, sten)););
598 }
599
600 Set::Matrix sig = Set::Matrix::Zero();
601
602 if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
603 sig = C(i, j, k).DW(gradu);
604 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
605 sig = C(i, j, k).DW(0.5 * (gradu + gradu.transpose()));
606 else if (C(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
607 sig = C(i, j, k).DW(gradu + Set::Matrix::Identity());
608
609 // = C(i,j,k)(gradu,m_homogeneous);
610
611 AMREX_D_PICK(dw(i, j, k, 0) = sig(0, 0);
612 ,
613 dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1);
614 dw(i, j, k, 2) = sig(1, 0); dw(i, j, k, 3) = sig(1, 1);
615 ,
616 dw(i, j, k, 0) = sig(0, 0); dw(i, j, k, 1) = sig(0, 1); dw(i, j, k, 2) = sig(0, 2);
617 dw(i, j, k, 3) = sig(1, 0); dw(i, j, k, 4) = sig(1, 1); dw(i, j, k, 5) = sig(1, 2);
618 dw(i, j, k, 6) = sig(2, 0); dw(i, j, k, 7) = sig(2, 1); dw(i, j, k, 8) = sig(2, 2););
619
620 });
621 }
622 }
623 }
624
625
626public:
627 int m_nriters = 1;
628 Set::Scalar m_nrtolerance = 0.0;
629 Operator::Elastic<T::sym>* m_elastic;
630 //BC::Operator::Elastic::Elastic *m_bc;
631
632 Set::Field<Set::Scalar>* m_psi = nullptr;
633
634public:
635 // These paramters control a standard Newton-Raphson solve.
636 //
637 // **Note**:
638 // This class inherits all of the linear solve paramters
639 // from its parent class, :ref:`Solver::Nonlocal::Linear`
640 static void Parse(Newton<T>& value, amrex::ParmParse& pp)
641 {
642 Linear::Parse(value, pp);
643
644 // Number of newton-raphson iterations.
645 pp_query("nriters", value.m_nriters);
646
647 // Tolerance to use for newton-raphson convergence
648 pp_query("nrtolerance", value.m_nrtolerance);
649 }
650
651};
652} // namespace Nonlocal
653} // namespace Solver
654
655
656#endif
pp_query
#define pp_query(...)
Definition ParmParse.H:106
TEST
#define TEST(x)
Definition Util.H:21
INFO
#define INFO
Definition Util.H:20
Operator::Elastic::SetModel
void SetModel(Set::Matrix4< AMREX_SPACEDIM, SYM > &a_model)
Definition Elastic.cpp:59
Operator::Elastic::GetBC
::BC::Operator::Elastic::Elastic & GetBC()
Definition Elastic.H:70
Operator::Elastic::SetPsi
void SetPsi(int amrlev, const amrex::MultiFab &a_psi)
Definition Elastic.cpp:126
Operator::Operator< Grid::Node >::Geom
const Geometry & Geom(int amr_lev, int mglev=0) const noexcept
Definition Operator.H:49
Set::Field::Copy
void Copy(int, amrex::MultiFab &, int, int) const
Definition Set.H:68
Set::Field::AddFrom
void AddFrom(int, amrex::MultiFab &, int, int) const
Definition Set.H:71
Set::Field::finest_level
int finest_level
Definition Set.H:67
Set::Field::Define
void Define(int a_levs, const amrex::Vector< amrex::BoxArray > &a_grids, const amrex::Vector< amrex::DistributionMapping > &a_dmap, int a_ncomp, int a_nghost)
Definition Set.H:52
Set::Field::NComp
int NComp() const
Definition Set.H:72
Set::Field::CopyFrom
void CopyFrom(int, amrex::MultiFab &, int, int) const
Definition Set.H:69
Solver::Nonlocal::Linear
Multigrid Linear solver for multicomponent, multi-level operators.
Definition Linear.H:19
Solver::Nonlocal::Linear::solve
Set::Scalar solve(amrex::Vector< std::unique_ptr< amrex::MultiFab > > &a_sol, amrex::Vector< std::unique_ptr< amrex::MultiFab > > &a_rhs, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Linear.H:132
Solver::Nonlocal::Linear::tol_rel
Set::Scalar tol_rel
Definition Linear.H:247
Solver::Nonlocal::Linear::Parse
static void Parse(Linear &value, amrex::ParmParse &pp)
Definition Linear.H:297
Solver::Nonlocal::Linear::Define
void Define(Operator::Operator< Grid::Node > &a_lp)
Definition Linear.H:36
Solver::Nonlocal::Linear::mlmg
amrex::MLMG * mlmg
Definition Linear.H:254
Solver::Nonlocal::Linear::tol_abs
Set::Scalar tol_abs
Definition Linear.H:248
Solver::Nonlocal::Linear::linop
Operator::Operator< Grid::Node > * linop
Definition Linear.H:253
Solver::Nonlocal::Newton::setPsi
void setPsi(Set::Field< Set::Scalar > &a_psi)
Definition Newton.H:46
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Scalar > &a_u_mf, const Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:435
Solver::Nonlocal::Newton::Newton
Newton(Operator::Elastic< T::sym > &a_op)
Definition Newton.H:20
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Scalar > &a_u_mf, const Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Newton.H:356
Solver::Nonlocal::Newton::m_nrtolerance
Set::Scalar m_nrtolerance
Definition Newton.H:628
Solver::Nonlocal::Newton::compResidual
void compResidual(Set::Field< Set::Vector > &a_res_mf, Set::Field< Set::Vector > &a_u_mf, Set::Field< Set::Vector > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:465
Solver::Nonlocal::Newton::Define
void Define(Operator::Elastic< T::sym > &a_op)
Definition Newton.H:30
Solver::Nonlocal::Newton::m_psi
Set::Field< Set::Scalar > * m_psi
Definition Newton.H:632
Solver::Nonlocal::Newton::W
void W(Set::Field< Set::Scalar > &a_w_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:525
Solver::Nonlocal::Newton::prepareForSolve
void prepareForSolve(const Set::Field< Set::Scalar > &a_u_mf, const Set::Field< Set::Scalar > &a_b_mf, Set::Field< Set::Scalar > &a_rhs_mf, Set::Field< Set::Matrix > &a_dw_mf, Set::Field< Set::Matrix4< AMREX_SPACEDIM, T::sym > > &a_ddw_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:52
Solver::Nonlocal::Newton::Parse
static void Parse(Newton< T > &value, amrex::ParmParse &pp)
Definition Newton.H:640
Solver::Nonlocal::Newton::m_elastic
Operator::Elastic< T::sym > * m_elastic
Definition Newton.H:629
Solver::Nonlocal::Newton::DW
void DW(Set::Field< Set::Scalar > &a_dw_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:568
Solver::Nonlocal::Newton::prepareForSolve
void prepareForSolve(const Set::Field< Set::Vector > &a_u_mf, const Set::Field< Set::Vector > &a_b_mf, Set::Field< Set::Scalar > &a_rhs_mf, Set::Field< Set::Matrix > &a_dw_mf, Set::Field< Set::Matrix4< AMREX_SPACEDIM, T::sym > > &a_ddw_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:144
Solver::Nonlocal::Newton::compLinearSolverResidual
void compLinearSolverResidual(Set::Field< Set::Vector > &a_res_mf, Set::Field< Set::Vector > &a_u_mf, Set::Field< Set::Vector > &a_b_mf)
Definition Newton.H:499
Solver::Nonlocal::Newton::setNRIters
void setNRIters(int a_nriters)
Definition Newton.H:44
Solver::Nonlocal::Newton::compResidual
void compResidual(Set::Field< Set::Scalar > &a_res_mf, Set::Field< Set::Scalar > &a_u_mf, Set::Field< Set::Scalar > &a_b_mf, Set::Field< T > &a_model_mf)
Definition Newton.H:442
Solver::Nonlocal::Newton::solve
Set::Scalar solve(const Set::Field< Set::Vector > &a_u_mf, const Set::Field< Set::Vector > &a_b_mf, Set::Field< T > &a_model_mf, Real a_tol_rel, Real a_tol_abs, const char *checkpoint_file=nullptr)
Definition Newton.H:275
Numeric::GetStencil
static AMREX_FORCE_INLINE std::array< StencilType, AMREX_SPACEDIM > GetStencil(const int i, const int j, const int k, const amrex::Box domain)
Definition Stencil.H:36
Numeric::CellGradientOnNode
AMREX_FORCE_INLINE Set::Vector CellGradientOnNode(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM])
Definition Stencil.H:690
Numeric::Gradient
AMREX_FORCE_INLINE Set::Vector Gradient(const amrex::Array4< const Set::Scalar > &f, const int &i, const int &j, const int &k, const int &m, const Set::Scalar dx[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition Stencil.H:672
Numeric::Divergence
AMREX_FORCE_INLINE Set::Vector Divergence(const amrex::Array4< const Set::Matrix > &dw, const int &i, const int &j, const int &k, const Set::Scalar DX[AMREX_SPACEDIM], std::array< StencilType, AMREX_SPACEDIM > &stencil=DefaultType)
Definition Stencil.H:581
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:20
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:23
Solver
A bunch of solvers.
Definition CG.H:6
Util::RealFillBoundary
AMREX_FORCE_INLINE void RealFillBoundary(amrex::FabArray< amrex::BaseFab< T > > &a_mf, const amrex::Geometry &)
Definition Util.H:322
Util::Assert
AMREX_FORCE_INLINE void Assert(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:70
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:141
Numeric::Interpolate::CellToNodeAverage
static AMREX_FORCE_INLINE T CellToNodeAverage(const amrex::Array4< const T > &f, const int &i, const int &j, const int &k, const int &m, std::array< StencilType, AMREX_SPACEDIM > stencil=DefaultType)
Definition Stencil.H:1293
Generated by doxygen 1.9.8