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 25 : Newton() {};
19 :
20 0 : Newton(Operator::Elastic<T::sym>& a_op)
21 0 : {
22 0 : this->Define(a_op);
23 0 : };
24 :
25 25 : ~Newton()
26 : {
27 25 : if (m_defined) Clear();
28 25 : }
29 :
30 426 : void Define(Operator::Elastic<T::sym>& a_op)
31 : {
32 426 : Linear::Define(a_op);
33 426 : m_elastic = dynamic_cast<Operator::Elastic<T::sym> *>(linop);
34 : //m_bc = &m_elastic->GetBC();
35 426 : }
36 426 : void Clear()
37 : {
38 426 : Linear::Clear();
39 426 : m_elastic = nullptr;
40 : //m_bc = nullptr;
41 426 : }
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 0 : 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 0 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
60 : {
61 0 : amrex::Box domain(linop->Geom(lev).Domain());
62 0 : domain.convert(amrex::IntVect::TheNodeVector());
63 0 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
64 0 : Set::Vector DX(linop->Geom(lev).CellSize());
65 0 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
66 : {
67 0 : amrex::Box bx = mfi.grownnodaltilebox();
68 0 : bx = bx & domain;
69 :
70 0 : amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
71 0 : amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
72 0 : amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
73 0 : 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 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
77 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
78 :
79 0 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
80 :
81 0 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
82 : {
83 0 : dw(i, j, k) = model(i, j, k).DW(gradu);
84 0 : ddw(i, j, k) = model(i, j, k).DDW(gradu);
85 : }
86 0 : else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::epsilon)
87 : {
88 0 : Set::Matrix eps = 0.5 * (gradu + gradu.transpose());
89 0 : dw(i, j, k) = model(i, j, k).DW(eps);
90 0 : ddw(i, j, k) = model(i, j, k).DDW(eps);
91 : }
92 0 : else if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
93 : {
94 0 : Set::Matrix F = gradu + Set::Matrix::Identity();
95 0 : dw(i, j, k) = model(i, j, k).DW(F);
96 0 : ddw(i, j, k) = model(i, j, k).DDW(F);
97 : }
98 : });
99 : }
100 :
101 0 : Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
102 0 : Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
103 : }
104 :
105 0 : m_elastic->SetModel(a_ddw_mf);
106 :
107 0 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
108 : {
109 0 : amrex::Box domain(linop->Geom(lev).Domain());
110 0 : domain.convert(amrex::IntVect::TheNodeVector());
111 0 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
112 0 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
113 0 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
114 : {
115 0 : amrex::Box bx = mfi.nodaltilebox();
116 0 : bx = bx & domain;
117 0 : amrex::Array4<const Set::Scalar> const& u = a_u_mf[lev]->array(mfi);
118 0 : amrex::Array4<const Set::Scalar> const& b = a_b_mf[lev]->array(mfi);
119 0 : amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
120 0 : amrex::Array4<Set::Scalar> const& rhs = a_rhs_mf[lev]->array(mfi);
121 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
122 : {
123 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
124 : // Do this if on the domain boundary
125 0 : if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
126 : {
127 0 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
128 0 : Set::Vector U(AMREX_D_DECL(u(i, j, k, 0), u(i, j, k, 1), u(i, j, k, 2)));
129 0 : Set::Vector ret = m_elastic->GetBC()(U, gradu, dw(i, j, k), i, j, k, bx);
130 0 : 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 0 : Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx);
135 0 : for (int p = 0; p < AMREX_SPACEDIM; p++) rhs(i, j, k, p) = b(i, j, k, p) - divdw(p);
136 : }
137 : });
138 : }
139 0 : Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
140 0 : Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
141 : }
142 0 : }
143 :
144 1326 : 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 2698 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
152 : {
153 1372 : amrex::Box domain(linop->Geom(lev).Domain());
154 1372 : domain.convert(amrex::IntVect::TheNodeVector());
155 1372 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
156 1372 : Set::Vector DX(linop->Geom(lev).CellSize());
157 2877 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
158 : {
159 1505 : amrex::Box bx = mfi.grownnodaltilebox();
160 1505 : bx = bx & domain;
161 :
162 1505 : amrex::Array4<const T> const& model = a_model_mf[lev]->array(mfi);
163 1505 : amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
164 1505 : amrex::Array4<Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
165 1505 : 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 617906 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
169 : {
170 209757 : auto sten = Numeric::GetStencil(i, j, k, bx);
171 :
172 209757 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
173 209757 : Set::Matrix kinvar;
174 419514 : if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::gradu)
175 138723 : 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 629271 : dw(i, j, k) = model(i, j, k).DW(kinvar);
182 629271 : ddw(i, j, k) = model(i, j, k).DDW(kinvar);
183 :
184 : });
185 : }
186 :
187 1372 : Util::RealFillBoundary(*a_dw_mf[lev], m_elastic->Geom(lev));
188 1372 : Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
189 : }
190 :
191 1326 : m_elastic->SetModel(a_ddw_mf);
192 1326 : if (m_psi)
193 58 : for (int i = 0; i <= a_b_mf.finest_level; i++)
194 44 : m_elastic->SetPsi(i, *(*m_psi)[i]);
195 :
196 :
197 2698 : for (int lev = 0; lev <= a_b_mf.finest_level; ++lev)
198 : {
199 1372 : amrex::Box domain(linop->Geom(lev).Domain());
200 1372 : domain.convert(amrex::IntVect::TheNodeVector());
201 1372 : const Set::Scalar* dx = linop->Geom(lev).CellSize();
202 1372 : const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
203 2877 : for (MFIter mfi(*a_model_mf[lev], false); mfi.isValid(); ++mfi)
204 : {
205 1505 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
206 1505 : amrex::Array4<const Set::Vector> const& u = a_u_mf[lev]->array(mfi);
207 1505 : amrex::Array4<const Set::Vector> const& b = a_b_mf[lev]->array(mfi);
208 1505 : amrex::Array4<const Set::Matrix> const& dw = a_dw_mf[lev]->array(mfi);
209 1505 : 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 1505 : if (m_psi == nullptr)
213 : {
214 122666 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
215 : {
216 121275 : auto sten = Numeric::GetStencil(i, j, k, bx);
217 : // Do this if on the domain boundary
218 121275 : if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
219 : {
220 46805 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
221 140415 : Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k), i, j, k, bx);
222 292825 : 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 74470 : Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
227 388550 : 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 114 : const amrex::Dim3 boxlo = amrex::lbound(bx), boxhi = amrex::ubound(bx);
237 114 : amrex::Array4<const Set::Scalar> const& psi = (*m_psi)[lev]->array(mfi);
238 88596 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
239 : {
240 88482 : auto sten = Numeric::GetStencil(i, j, k, bx);
241 : // Do this if on the domain boundary
242 88482 : if (AMREX_D_TERM(i == lo.x || i == hi.x, || j == lo.y || j == hi.y, || k == lo.z || k == hi.z))
243 : {
244 2560 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, dx, sten);
245 2560 : Set::Scalar psiavg = 1.0;
246 2560 : if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
247 : {
248 0 : psiavg = Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
249 : }
250 7680 : Set::Vector ret = m_elastic->GetBC()(u(i, j, k), gradu, dw(i, j, k) * psiavg, i, j, k, bx);
251 12800 : 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 85922 : Set::Vector divdw = Numeric::Divergence(dw, i, j, k, dx, sten);
256 85922 : if (AMREX_D_TERM(i > boxlo.x && i<boxhi.x, && j>boxlo.y && j<boxhi.y, && k>boxlo.z && k < boxhi.z))
257 : {
258 153860 : divdw *= Numeric::Interpolate::CellToNodeAverage(psi, i, j, k, 0);
259 76930 : Set::Vector gradpsi = Numeric::CellGradientOnNode(psi, i, j, k, 0, dx);
260 153860 : divdw += dw(i, j, k) * gradpsi;
261 : }
262 429610 : 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 1372 : Util::RealFillBoundary(*a_ddw_mf[lev], m_elastic->Geom(lev));
269 1372 : Util::RealFillBoundary(*a_rhs_mf[lev], m_elastic->Geom(lev));
270 : }
271 1326 : }
272 :
273 :
274 : public:
275 426 : 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 852 : Set::Field<Set::Scalar> dsol_mf, rhs_mf;
281 852 : Set::Field<Set::Matrix> dw_mf;
282 852 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
283 :
284 426 : dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
285 426 : dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
286 426 : ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
287 426 : rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
288 898 : for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
289 : {
290 944 : dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
291 472 : a_u_mf[lev]->DistributionMap(),
292 : a_u_mf.NComp(),
293 472 : a_u_mf[lev]->nGrow());
294 944 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
295 472 : a_b_mf[lev]->DistributionMap(),
296 : 1,
297 472 : a_b_mf[lev]->nGrow());
298 944 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
299 472 : a_b_mf[lev]->DistributionMap(),
300 : 1,
301 472 : a_b_mf[lev]->nGrow());
302 944 : rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
303 472 : a_b_mf[lev]->DistributionMap(),
304 : a_b_mf.NComp(),
305 472 : a_b_mf[lev]->nGrow());
306 :
307 472 : dsol_mf[lev]->setVal(0.0);
308 472 : dw_mf[lev]->setVal(Set::Matrix::Zero());
309 472 : ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
310 :
311 472 : 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 1750 : for (int nriter = 0; nriter < m_nriters; nriter++)
316 : {
317 1326 : if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
318 :
319 1326 : prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
320 :
321 :
322 1326 : if (nriter == m_nriters) break;
323 :
324 1326 : Solver::Nonlocal::Linear::solve(dsol_mf, rhs_mf, a_tol_rel, a_tol_abs, checkpoint_file);
325 :
326 1326 : Set::Scalar cornorm = 0, solnorm = 0;
327 2698 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
328 : {
329 4416 : for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
330 : {
331 3044 : Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
332 3044 : 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 1326 : if (solnorm == 0) relnorm = cornorm;
341 0 : else relnorm = cornorm / solnorm;
342 1326 : if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
343 :
344 2698 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
345 1372 : 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 1326 : if (relnorm < m_nrtolerance)
349 2 : return relnorm;
350 :
351 : }
352 :
353 424 : return 0.0;
354 : }
355 :
356 0 : 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 0 : Set::Field<Set::Scalar> dsol_mf, rhs_mf;
362 0 : Set::Field<Set::Matrix> dw_mf;
363 0 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
364 :
365 0 : dsol_mf.resize(a_u_mf.finest_level + 1); dsol_mf.finest_level = a_u_mf.finest_level;
366 0 : dw_mf.resize(a_u_mf.finest_level + 1); dw_mf.finest_level = a_u_mf.finest_level;
367 0 : ddw_mf.resize(a_u_mf.finest_level + 1); ddw_mf.finest_level = a_u_mf.finest_level;
368 0 : rhs_mf.resize(a_u_mf.finest_level + 1); rhs_mf.finest_level = a_u_mf.finest_level;
369 0 : for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
370 : {
371 0 : dsol_mf.Define(lev, a_u_mf[lev]->boxArray(),
372 0 : a_u_mf[lev]->DistributionMap(),
373 0 : a_u_mf[lev]->nComp(),
374 0 : a_u_mf[lev]->nGrow());
375 0 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
376 0 : a_b_mf[lev]->DistributionMap(),
377 : 1,
378 0 : a_b_mf[lev]->nGrow());
379 0 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
380 0 : a_b_mf[lev]->DistributionMap(),
381 : 1,
382 0 : a_b_mf[lev]->nGrow());
383 0 : rhs_mf.Define(lev, a_b_mf[lev]->boxArray(),
384 0 : a_b_mf[lev]->DistributionMap(),
385 0 : a_b_mf[lev]->nComp(),
386 0 : a_b_mf[lev]->nGrow());
387 :
388 0 : dsol_mf[lev]->setVal(0.0);
389 0 : dw_mf[lev]->setVal(Set::Matrix::Zero());
390 0 : ddw_mf[lev]->setVal(Set::Matrix4<AMREX_SPACEDIM, T::sym>::Zero());
391 :
392 0 : amrex::MultiFab::Copy(*rhs_mf[lev], *a_b_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
393 : }
394 :
395 0 : for (int nriter = 0; nriter < m_nriters; nriter++)
396 : {
397 0 : if (verbose > 0 && nriter < m_nriters) Util::Message(INFO, "Newton Iteration ", nriter + 1, " of ", m_nriters);
398 :
399 0 : prepareForSolve(a_u_mf, a_b_mf, rhs_mf, dw_mf, ddw_mf, a_model_mf);
400 :
401 :
402 0 : if (nriter == m_nriters) break;
403 :
404 0 : 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 0 : Set::Scalar cornorm = 0, solnorm = 0;
408 0 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
409 : {
410 0 : for (int comp = 0; comp < AMREX_SPACEDIM; comp++)
411 : {
412 0 : Set::Scalar tmpcornorm = dsol_mf[lev]->norm0(comp, 0);
413 0 : if (tmpcornorm > cornorm) cornorm = tmpcornorm;
414 :
415 0 : Set::Scalar tmpsolnorm = a_u_mf[lev]->norm0(comp, 0);
416 0 : if (tmpsolnorm > solnorm) solnorm = tmpsolnorm;
417 : }
418 :
419 : }
420 : Set::Scalar relnorm;
421 0 : if (solnorm == 0) relnorm = cornorm;
422 0 : else relnorm = cornorm / solnorm;
423 0 : if (verbose > 0) Util::Message(INFO, "NR iteration ", nriter + 1, ", relative norm(ddisp) = ", relnorm);
424 :
425 0 : for (int lev = 0; lev < dsol_mf.size(); ++lev)
426 0 : amrex::MultiFab::Add(*a_u_mf[lev], *dsol_mf[lev], 0, 0, AMREX_SPACEDIM, 2);
427 :
428 0 : if (relnorm < m_nrtolerance)
429 0 : return relnorm;
430 :
431 : }
432 :
433 0 : 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 0 : 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 0 : Set::Field<Set::Matrix> dw_mf;
448 0 : Set::Field<Set::Matrix4<AMREX_SPACEDIM, T::sym>> ddw_mf;
449 0 : dw_mf.resize(a_u_mf.size());
450 0 : ddw_mf.resize(a_u_mf.size());
451 0 : for (int lev = 0; lev < a_u_mf.size(); lev++)
452 : {
453 0 : dw_mf.Define(lev, a_b_mf[lev]->boxArray(),
454 0 : a_b_mf[lev]->DistributionMap(),
455 0 : 1, a_b_mf[lev]->nGrow());
456 0 : ddw_mf.Define(lev, a_b_mf[lev]->boxArray(),
457 0 : a_b_mf[lev]->DistributionMap(),
458 0 : 1, a_b_mf[lev]->nGrow());
459 0 : dw_mf[lev]->setVal(Set::Matrix::Zero());
460 : }
461 :
462 0 : prepareForSolve(a_u_mf, a_b_mf, a_res_mf, dw_mf, ddw_mf, a_model_mf);
463 0 : }
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 4 : 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 8 : Util::Assert(INFO, TEST(a_res_mf.finest_level == a_u_mf.finest_level));
504 8 : Util::Assert(INFO, TEST(a_u_mf.finest_level == a_b_mf.finest_level));
505 8 : 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 18 : for (int lev = 0; lev <= a_u_mf.finest_level; lev++)
507 : {
508 14 : res_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
509 14 : sol_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
510 14 : rhs_mf.Define(lev, a_b_mf[lev]->boxArray(), a_b_mf[lev]->DistributionMap(), AMREX_SPACEDIM, a_b_mf[lev]->nGrow());
511 :
512 14 : a_u_mf.Copy(lev, *sol_mf[lev], 0, 2);
513 14 : a_b_mf.Copy(lev, *rhs_mf[lev], 0, 2);
514 : }
515 :
516 4 : mlmg->compResidual(amrex::GetVecOfPtrs(res_mf), amrex::GetVecOfPtrs(sol_mf), amrex::GetVecOfConstPtrs(rhs_mf));
517 :
518 18 : for (int lev = 0; lev <= a_res_mf.finest_level; ++lev)
519 : {
520 14 : a_res_mf.CopyFrom(lev, *res_mf[lev], 0, 2);
521 : }
522 4 : }
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 :
626 : public:
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 :
634 : public:
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 19 : static void Parse(Newton<T>& value, amrex::ParmParse& pp)
641 : {
642 19 : Linear::Parse(value, pp);
643 :
644 : // Number of newton-raphson iterations.
645 19 : pp_query("nriters", value.m_nriters);
646 :
647 : // Tolerance to use for newton-raphson convergence
648 19 : pp_query("nrtolerance", value.m_nrtolerance);
649 19 : }
650 :
651 : };
652 : } // namespace Nonlocal
653 : } // namespace Solver
654 :
655 :
656 : #endif
|
