Alamo
Hydro.cpp
Go to the documentation of this file.
1
2#include "Hydro.H"
3#include "IO/ParmParse.H"
4#include "BC/Constant.H"
5#include "BC/Expression.H"
6#include "Numeric/Stencil.H"
7#include "IC/Constant.H"
8#include "IC/Laminate.H"
9#include "IC/Expression.H"
10#include "IC/BMP.H"
11#include "IC/PNG.H"
12#include "Solver/Local/Riemann/Roe.H"
13#include "Solver/Local/Riemann/HLLE.H"
14#include "Solver/Local/Riemann/HLLC.H"
15#if AMREX_SPACEDIM == 2
16
17namespace Integrator
18{
19
20Hydro::Hydro(IO::ParmParse& pp) : Hydro()
21{
22 pp_queryclass(*this);
23}
24
25void
26Hydro::Parse(Hydro& value, IO::ParmParse& pp)
27{
28 BL_PROFILE("Integrator::Hydro::Hydro()");
29 {
30 // pp.query_default("r_refinement_criterion", value.r_refinement_criterion , 0.01);
31 // energy-based refinement
32 // pp.query_default("e_refinement_criterion", value.e_refinement_criterion , 0.01);
33 // momentum-based refinement
34 // pp.query_default("m_refinement_criterion", value.m_refinement_criterion , 0.01);
35
36 std::string scheme_str;
37 // time integration scheme to use
38 pp.query_validate("scheme",scheme_str, {"forwardeuler","ssprk3","rk4"});
39 if (scheme_str == "forwardeuler") value.scheme = IntegrationScheme::ForwardEuler;
40 else if (scheme_str == "ssprk3") value.scheme = IntegrationScheme::SSPRK3;
41 else if (scheme_str == "rk4") value.scheme = IntegrationScheme::RK4;
42
43
44 // eta-based refinement
45 pp.query_default("eta_refinement_criterion", value.eta_refinement_criterion , 0.01);
46 // vorticity-based refinement
47 pp.query_default("omega_refinement_criterion", value.omega_refinement_criterion, 0.01);
48 // velocity gradient-based refinement
49 pp.query_default("gradu_refinement_criterion", value.gradu_refinement_criterion, 0.01);
50 // pressure-based refinement
51 pp.query_default("p_refinement_criterion", value.p_refinement_criterion, 1e100);
52 // density-based refinement
53 pp.query_default("rho_refinement_criterion", value.rho_refinement_criterion, 1e100);
54
55 pp_query_required("gamma", value.gamma); // gamma for gamma law
56 pp_query_required("cfl", value.cfl); // cfl condition
57 pp_query_default("cfl_v", value.cfl_v,1E100); // cfl condition
58 pp_query_required("mu", value.mu); // linear viscosity coefficient
59 pp_forbid("Lfactor","replaced with mu");
60 //pp_query_default("Lfactor", value.Lfactor,1.0); // (to be removed) test factor for viscous source
61 pp_forbid("Pfactor","replaced with mu");
62 //pp_query_default("Pfactor", value.Pfactor,1.0); // (to be removed) test factor for viscous source
63 pp_query_default("pref", value.pref,1.0); // reference pressure for Roe solver
64
65 pp_forbid("rho.bc","--> density.bc");
66 pp_forbid("p.bc","--> pressure.bc");
67 pp_forbid("v.bc", "--> velocity.bc");
68 pp_forbid("pressure.bc","--> energy.bc");
69 pp_forbid("velocity.bc","--> momentum.bc");
70
71 // Boundary condition for density
72 pp.select_default<BC::Constant,BC::Expression>("density.bc",value.density_bc,1);
73 // Boundary condition for energy
74 pp.select_default<BC::Constant,BC::Expression>("energy.bc",value.energy_bc,1);
75 // Boundary condition for momentum
76 pp.select_default<BC::Constant,BC::Expression>("momentum.bc",value.momentum_bc,2);
77 // Boundary condition for phase field order parameter
78 pp.select_default<BC::Constant,BC::Expression>("pf.eta.bc",value.eta_bc,1);
79
80 pp_query_default("small",value.small,1E-8); // small regularization value
81 pp_query_default("cutoff",value.cutoff,-1E100); // cutoff value
82 pp_query_default("lagrange",value.lagrange,0.0); // lagrange no-penetration factor
83
84 pp_forbid("roefix","--> solver.roe.entropy_fix"); // Roe solver entropy fix
85
86 }
87 // Register FabFields:
88 {
89 int nghost = 1;
90
91 value.RegisterNewFab(value.eta_mf, value.eta_bc, 1, nghost, "eta", true, true);
92 value.RegisterNewFab(value.eta_old_mf, value.eta_bc, 1, nghost, "eta_old", true, true);
93 value.RegisterNewFab(value.etadot_mf, value.eta_bc, 1, nghost, "etadot", true, false);
94
95 value.RegisterNewFab(value.density_mf, value.density_bc, 1, nghost, "density", true , true);
96 value.RegisterNewFab(value.density_old_mf, value.density_bc, 1, nghost, "density_old", false, true);
97
98 value.RegisterNewFab(value.energy_mf, value.energy_bc, 1, nghost, "energy", true ,true);
99 value.RegisterNewFab(value.energy_old_mf, value.energy_bc, 1, nghost, "energy_old" , false, true);
100
101 value.RegisterNewFab(value.momentum_mf, value.momentum_bc, 2, nghost, "momentum", true ,true, {"x","y"});
102 value.RegisterNewFab(value.momentum_old_mf, value.momentum_bc, 2, nghost, "momentum_old", false, true);
103
104 value.RegisterNewFab(value.pressure_mf, &value.bc_nothing, 1, nghost, "pressure", true, false);
105 value.RegisterNewFab(value.velocity_mf, &value.bc_nothing, 2, nghost, "velocity", true, false,{"x","y"});
106 value.RegisterNewFab(value.vorticity_mf, &value.bc_nothing, 1, nghost, "vorticity", true, false);
107
108 value.RegisterNewFab(value.m0_mf, &value.bc_nothing, 1, 0, "m0", true, false);
109 value.RegisterNewFab(value.u0_mf, &value.bc_nothing, 2, 0, "u0", true, false, {"x","y"});
110 value.RegisterNewFab(value.q_mf, &value.bc_nothing, 2, 0, "q", true, false, {"x","y"});
111
112 value.RegisterNewFab(value.solid.momentum_mf, &value.neumann_bc_D, 2, nghost, "solid.momentum", true, false, {"x","y"});
113 value.RegisterNewFab(value.solid.density_mf, &value.neumann_bc_1, 1, nghost, "solid.density", true, false);
114 value.RegisterNewFab(value.solid.energy_mf, &value.neumann_bc_1, 1, nghost, "solid.energy", true, false);
115
116 value.RegisterNewFab(value.Source_mf, &value.bc_nothing, 4, 0, "Source", true, false);
117 }
118
119 pp_forbid("Velocity.ic.type", "--> velocity.ic.type");
120 pp_forbid("Pressure.ic", "--> pressure.ic");
121 pp_forbid("SolidMomentum.ic", "--> solid.momentum.ic");
122 pp_forbid("SolidDensity.ic.type", "--> solid.density.ic.type");
123 pp_forbid("SolidEnergy.ic.type", "--> solid.energy.ic.type");
124 pp_forbid("Density.ic.type", "--> density.ic.type");
125 pp_forbid("rho_injected.ic.type","no longer using rho_injected use m0 instead");
126 pp.forbid("mdot.ic.type", "replace mdot with u0");
127
128
129 // ORDER PARAMETER
130
131 // eta initial condition
132 pp.select_default<IC::Constant,IC::Laminate,IC::Expression,IC::BMP,IC::PNG>("eta.ic",value.eta_ic,value.geom);
133
134 // PRIMITIVE FIELD INITIAL CONDITIONS
135
136 // velocity initial condition
137 pp.select_default<IC::Constant,IC::Expression>("velocity.ic",value.velocity_ic,value.geom);
138 // solid pressure initial condition
139 pp.select_default<IC::Constant,IC::Expression>("pressure.ic",value.pressure_ic,value.geom);
140 // density initial condition type
141 pp.select_default<IC::Constant,IC::Expression>("density.ic",value.density_ic,value.geom);
142
143
144 // SOLID FIELDS
145
146 // solid momentum initial condition
147 pp.select_default<IC::Constant,IC::Expression>("solid.momentum.ic",value.solid.momentum_ic,value.geom);
148 // solid density initial condition
149 pp.select_default<IC::Constant,IC::Expression>("solid.density.ic",value.solid.density_ic,value.geom);
150 // solid energy initial condition
151 pp.select_default<IC::Constant,IC::Expression>("solid.energy.ic",value.solid.energy_ic,value.geom);
152
153
154 // DIFFUSE BOUNDARY SOURCES
155
156 // diffuse boundary prescribed mass flux
157 pp.select_default<IC::Constant,IC::Expression>("m0.ic",value.ic_m0,value.geom);
158 // diffuse boundary prescribed velocity
159 pp.select_default<IC::Constant,IC::Expression>("u0.ic",value.ic_u0,value.geom);
160 // diffuse boundary prescribed heat flux
161 pp.select_default<IC::Constant,IC::Expression>("q.ic",value.ic_q,value.geom);
162
163 // Riemann solver
164 pp.select_default< Solver::Local::Riemann::Roe,
165 Solver::Local::Riemann::HLLE,
166 Solver::Local::Riemann::HLLC>("solver",value.riemannsolver);
167
168
169
170 bool allow_unused;
171 // Set this to true to allow unused inputs without error.
172 // (Not recommended.)
173 pp.query_default("allow_unused",allow_unused,false);
174 if (!allow_unused && pp.AnyUnusedInputs(true, false))
175 {
176 Util::Warning(INFO,"The following inputs were specified but not used:");
177 pp.AllUnusedInputs();
178 Util::Exception(INFO,"Aborting. Specify 'allow_unused=True` to ignore this error.");
179 }
180}
181
182
183void Hydro::Initialize(int lev)
184{
185 BL_PROFILE("Integrator::Hydro::Initialize");
186
187 eta_ic ->Initialize(lev, eta_mf, 0.0);
188 eta_ic ->Initialize(lev, eta_old_mf, 0.0);
189 etadot_mf[lev] ->setVal(0.0);
190
191 //flux_mf[lev] ->setVal(0.0);
192
193 velocity_ic ->Initialize(lev, velocity_mf, 0.0);
194 pressure_ic ->Initialize(lev, pressure_mf, 0.0);
195 density_ic ->Initialize(lev, density_mf, 0.0);
196
197 density_ic ->Initialize(lev, density_old_mf, 0.0);
198
199
200 solid.density_ic ->Initialize(lev, solid.density_mf, 0.0);
201 solid.momentum_ic->Initialize(lev, solid.momentum_mf, 0.0);
202 solid.energy_ic ->Initialize(lev, solid.energy_mf, 0.0);
203
204 ic_m0 ->Initialize(lev, m0_mf, 0.0);
205 ic_u0 ->Initialize(lev, u0_mf, 0.0);
206 ic_q ->Initialize(lev, q_mf, 0.0);
207
208 Source_mf[lev] ->setVal(0.0);
209
210 Mix(lev);
211}
212
213void Hydro::Mix(int lev)
214{
215 for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
216 {
217 const amrex::Box& bx = mfi.growntilebox();
218
219 Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
220
221 Set::Patch<const Set::Scalar> v = velocity_mf.Patch(lev,mfi);
222 Set::Patch<const Set::Scalar> p = pressure_mf.Patch(lev,mfi);
223 Set::Patch<Set::Scalar> rho = density_mf.Patch(lev,mfi);
224 Set::Patch<Set::Scalar> rho_old = density_old_mf.Patch(lev,mfi);
225 Set::Patch<Set::Scalar> M = momentum_mf.Patch(lev,mfi);
226 Set::Patch<Set::Scalar> M_old = momentum_old_mf.Patch(lev,mfi);
227 Set::Patch<Set::Scalar> E = energy_mf.Patch(lev,mfi);
228 Set::Patch<Set::Scalar> E_old = energy_old_mf.Patch(lev,mfi);
229 Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
230 Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
231 Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
232
233
234 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
235 {
236 rho(i, j, k) = eta(i, j, k) * rho(i, j, k) + (1.0 - eta(i, j, k)) * rho_solid(i, j, k);
237 rho_old(i, j, k) = rho(i, j, k);
238
239 M(i, j, k, 0) = (rho(i, j, k)*v(i, j, k, 0))*eta(i, j, k) + M_solid(i, j, k, 0)*(1.0-eta(i, j, k));
240 M(i, j, k, 1) = (rho(i, j, k)*v(i, j, k, 1))*eta(i, j, k) + M_solid(i, j, k, 1)*(1.0-eta(i, j, k));
241 M_old(i, j, k, 0) = M(i, j, k, 0);
242 M_old(i, j, k, 1) = M(i, j, k, 1);
243
244 E(i, j, k) =
245 (0.5 * (v(i, j, k, 0) * v(i, j, k, 0) + v(i, j, k, 1) * v(i, j, k, 1)) * rho(i, j, k) + p(i, j, k) / (gamma - 1.0)) * eta(i, j, k)
246 +
247 E_solid(i, j, k) * (1.0 - eta(i, j, k));
248 E_old(i, j, k) = E(i, j, k);
249 });
250 }
251 c_max = 0.0;
252 vx_max = 0.0;
253 vy_max = 0.0;
254}
255
256void Hydro::UpdateEta(int lev, Set::Scalar time)
257{
258 eta_ic->Initialize(lev, eta_mf, time);
259}
260
261void Hydro::TimeStepBegin(Set::Scalar, int /*iter*/)
262{
263
264}
265
266void Hydro::TimeStepComplete(Set::Scalar, int lev)
267{
268 if (dynamictimestep.on)
269 Integrator::DynamicTimestep_Update();
270 return;
271
272 const Set::Scalar* DX = geom[lev].CellSize();
273
274 amrex::ParallelDescriptor::ReduceRealMax(c_max);
275 amrex::ParallelDescriptor::ReduceRealMax(vx_max);
276 amrex::ParallelDescriptor::ReduceRealMax(vy_max);
277
278 Set::Scalar new_timestep = cfl / ((c_max + vx_max) / DX[0] + (c_max + vy_max) / DX[1]);
279
280 Util::Assert(INFO, TEST(AMREX_SPACEDIM == 2));
281
282 SetTimestep(new_timestep);
283}
284
285void Hydro::Advance(int lev, Set::Scalar time, Set::Scalar dt)
286{
287
288 std::swap(eta_old_mf, eta_mf);
289 std::swap(density_old_mf[lev], density_mf[lev]);
290 std::swap(momentum_old_mf[lev], momentum_mf[lev]);
291 std::swap(energy_old_mf[lev], energy_mf[lev]);
292 Set::Scalar dt_max = std::numeric_limits<Set::Scalar>::max();
293
294 UpdateEta(lev, time);
295
296 for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
297 {
298 const amrex::Box& bx = mfi.growntilebox();
299 amrex::Array4<const Set::Scalar> const& eta_new = (*eta_mf[lev]).array(mfi);
300 amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
301 amrex::Array4<Set::Scalar> const& etadot = (*etadot_mf[lev]).array(mfi);
302 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
303 {
304 etadot(i, j, k) = (eta_new(i, j, k) - eta(i, j, k)) / dt;
305 });
306 }
307
308
309 if (scheme == IntegrationScheme::ForwardEuler) // forward euler
310 {
311
312 RHS(lev, time,
313 *density_mf[lev], *momentum_mf[lev], *energy_mf[lev],
314 *density_old_mf[lev], *momentum_old_mf[lev], *energy_old_mf[lev]);
315
316 for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
317 {
318 const amrex::Box& bx = mfi.validbox();
319
320 Set::Patch<const Set::Scalar> rho_rhs = density_mf.Patch(lev,mfi);
321 Set::Patch<const Set::Scalar> E_rhs = energy_mf.Patch(lev,mfi);
322 Set::Patch<const Set::Scalar> M_rhs = momentum_mf.Patch(lev,mfi);
323
324 Set::Patch<const Set::Scalar> rho_old = density_old_mf.Patch(lev,mfi);
325 Set::Patch<const Set::Scalar> E_old = energy_old_mf.Patch(lev,mfi);
326 Set::Patch<const Set::Scalar> M_old = momentum_old_mf.Patch(lev,mfi);
327
328 Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
329 Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
330 Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
331
332
333 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
334 {
335 rho_new(i, j, k) = rho_old(i, j, k) + dt * rho_rhs(i,j,k);
336 M_new(i,j,k,0) = M_old(i,j,k,0) + dt * M_rhs(i,j,k,0);
337 M_new(i,j,k,1) = M_old(i,j,k,1) + dt * M_rhs(i,j,k,1);
338 E_new(i,j,k) = E_old(i,j,k) + dt * E_rhs(i,j,k);
339 });
340 }
341 }
342
343
344 else if (scheme == IntegrationScheme::SSPRK3)
345 {
346 // Butcher Tableau
347 // |
348 // 1 | 1
349 // 1/2 | 1/4 1/4
350 // ---------------------
351 // | 1/6 1/6 2/3
352
353 Set::Scalar
354 /* */
355 /* */ c2 = 1.0 , a21 = 1.0,
356 /* */ c3 = 0.5, a31 = 0.25, a32 = 0.25,
357 /* --------------------------------------------- */
358 /* */ b1 = 1./6, b2 = 1./6., b3 = 2./3.;
359
360
361 const amrex::BoxArray &ba = density_mf[lev]->boxArray();
362 const amrex::DistributionMapping &dm = density_mf[lev]->DistributionMap();
363 const int ng = density_mf[lev]->nGrow();
364
365 // handles to old solution
366 const amrex::MultiFab &density_old = *density_old_mf[lev];
367 const amrex::MultiFab &momentum_old = *momentum_old_mf[lev];
368 const amrex::MultiFab &energy_old = *energy_old_mf[lev];
369
370
371 // temporary storage
372 amrex::MultiFab density_k1(ba,dm,1,0), momentum_k1(ba,dm,2,0), energy_k1(ba,dm,1,0);
373 amrex::MultiFab density_k2(ba,dm,1,0), momentum_k2(ba,dm,2,0), energy_k2(ba,dm,1,0);
374 amrex::MultiFab density_k3(ba,dm,1,0), momentum_k3(ba,dm,2,0), energy_k3(ba,dm,1,0);
375
376 // buffer to hold combs of k1
377 amrex::MultiFab density_temp(ba,dm,1,ng), momentum_temp(ba,dm,2,ng), energy_temp(ba,dm,1,ng);
378
379 // fill the ghost cells from the _old fields, which were updated from the coarse patch.
380 density_temp.ParallelCopyToGhost(*density_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
381 momentum_temp.ParallelCopyToGhost(*momentum_old_mf[lev],0,0,2,amrex::IntVect(1),amrex::IntVect(1));
382 energy_temp.ParallelCopyToGhost(*energy_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
383
384 // handles to new solution
385 amrex::MultiFab &density_new = *density_mf[lev];
386 amrex::MultiFab &momentum_new = *momentum_mf[lev];
387 amrex::MultiFab &energy_new = *energy_mf[lev];
388
389
390 //
391 // Calculate K1
392 //
393 // k1 = RHS(t, yold)
394
395 RHS(lev,time,
396 density_k1,momentum_k1,energy_k1,
397 *density_old_mf[lev], *momentum_old_mf[lev], *energy_old_mf[lev]);
398
399
400 //
401 // Calculate K2
402 //
403 // ytemp = yold + dt*( a21*k1 )
404 amrex::MultiFab::LinComb(density_temp, 1.0, density_old, 0, dt*a21, density_k1, 0, 0, 1, 0);
405 amrex::MultiFab::LinComb(momentum_temp, 1.0, momentum_old, 0, dt*a21, momentum_k1, 0, 0, 2, 0);
406 amrex::MultiFab::LinComb(energy_temp, 1.0, energy_old, 0, dt*a21, energy_k1, 0, 0, 1, 0);
407 // fill boundary
408 density_bc ->FillBoundary(density_temp, 0, 1, time, 0); density_temp.FillBoundary(true);
409 momentum_bc->FillBoundary(momentum_temp, 0, 2, time, 0); momentum_temp.FillBoundary(true);
410 energy_bc ->FillBoundary(energy_temp, 0, 1, time, 0); energy_temp.FillBoundary(true);
411 // k2 = RHS(t + c2*dt, ytemp)
412 RHS(lev,time + c2*dt,
413 density_k2, momentum_k2, energy_k2,
414 density_temp, momentum_temp, energy_temp);
415
416 //
417 // Calculate K3
418 //
419 // ytemp = yold + dt*( a31*k1 + a32*k2 )
420 //
421 // 1. ytemp = yold + dt*a31*k1
422 amrex::MultiFab::LinComb(density_temp, 1.0, density_old, 0, dt*a31, density_k1, 0, 0, 1, 0);
423 amrex::MultiFab::LinComb(momentum_temp, 1.0, momentum_old, 0, dt*a31, momentum_k1, 0, 0, 2, 0);
424 amrex::MultiFab::LinComb(energy_temp, 1.0, energy_old, 0, dt*a31, energy_k1, 0, 0, 1, 0);
425 // 2. ytemp += dt*a32*k2
426 amrex::MultiFab::Saxpy(density_temp, dt*a32, density_k2, 0, 0, 1, 0);
427 amrex::MultiFab::Saxpy(momentum_temp, dt*a32, momentum_k2, 0, 0, 2, 0);
428 amrex::MultiFab::Saxpy(energy_temp, dt*a32, energy_k2, 0, 0, 1, 0);
429 // 3. fill boundary
430 density_bc ->FillBoundary(density_temp, 0, 1, time+c2*dt, 0); density_temp.FillBoundary(true);
431 momentum_bc->FillBoundary(momentum_temp, 0, 2, time+c2*dt, 0); momentum_temp.FillBoundary(true);
432 energy_bc ->FillBoundary(energy_temp, 0, 1, time+c2*dt, 0); energy_temp.FillBoundary(true);
433 // 4. k3 = RHS(t + c3*dt, ytemp)
434 RHS(lev,time + c3*dt,
435 density_k3, momentum_k3, energy_k3,
436 density_temp, momentum_temp, energy_temp);
437
438 //
439 // Assemble to get ynew
440 //
441 // ynew = yold + dt*(b1*k1 + b2*k2 + b3*k3)
442 //
443 // 1. ynew = yold + dt*b1*k1
444 amrex::MultiFab::LinComb(density_new, 1.0, density_old, 0, dt*b1, density_k1, 0, 0, 1, 0);
445 amrex::MultiFab::LinComb(momentum_new, 1.0, momentum_old, 0, dt*b1, momentum_k1, 0, 0, 2, 0);
446 amrex::MultiFab::LinComb(energy_new, 1.0, energy_old, 0, dt*b1, energy_k1, 0, 0, 1, 0);
447 // 2. ynew += dt*b2*k2
448 amrex::MultiFab::Saxpy(density_new, dt*b2, density_k2, 0, 0, 1, 0);
449 amrex::MultiFab::Saxpy(momentum_new, dt*b2, momentum_k2, 0, 0, 2, 0);
450 amrex::MultiFab::Saxpy(energy_new, dt*b2, energy_k2, 0, 0, 1, 0);
451 // 2. ynew += dt*b3*k3
452 amrex::MultiFab::Saxpy(density_new, dt*b3, density_k3, 0, 0, 1, 0);
453 amrex::MultiFab::Saxpy(momentum_new, dt*b3, momentum_k3, 0, 0, 2, 0);
454 amrex::MultiFab::Saxpy(energy_new, dt*b3, energy_k3, 0, 0, 1, 0);
455
456 }
457
458
459
460 else if (scheme == IntegrationScheme::RK4)
461 {
462 //
463 // RK4 time integration scheme:
464 //
465
466 const amrex::BoxArray &ba = density_mf[lev]->boxArray();
467 const amrex::DistributionMapping &dm = density_mf[lev]->DistributionMap();
468 const int ng = density_mf[lev]->nGrow();
469
470 // handles to old solution
471 const amrex::MultiFab &density_old = *density_old_mf[lev];
472 const amrex::MultiFab &momentum_old = *momentum_old_mf[lev];
473 const amrex::MultiFab &energy_old = *energy_old_mf[lev];
474
475 // runge kutta stages
476 amrex::MultiFab density_k1(ba,dm,1,0), momentum_k1(ba,dm,2,0), energy_k1(ba,dm,1,0);
477 amrex::MultiFab density_k2(ba,dm,1,0), momentum_k2(ba,dm,2,0), energy_k2(ba,dm,1,0);
478 amrex::MultiFab density_k3(ba,dm,1,0), momentum_k3(ba,dm,2,0), energy_k3(ba,dm,1,0);
479 amrex::MultiFab density_k4(ba,dm,1,0), momentum_k4(ba,dm,2,0), energy_k4(ba,dm,1,0);
480
481 // temporary storage
482 amrex::MultiFab density_st(ba,dm,1,ng), momentum_st(ba,dm,2,ng), energy_st(ba,dm,1,ng);
483
484 // fill the ghost cells from the _old fields, which were updated from the coarse patch.
485 density_st.ParallelCopyToGhost(*density_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
486 momentum_st.ParallelCopyToGhost(*momentum_old_mf[lev],0,0,2,amrex::IntVect(1),amrex::IntVect(1));
487 energy_st.ParallelCopyToGhost(*energy_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
488
489
490 // handles to new solution
491 amrex::MultiFab &density_new = *density_mf[lev];
492 amrex::MultiFab &momentum_new = *momentum_mf[lev];
493 amrex::MultiFab &energy_new = *energy_mf[lev];
494
495
496 //
497 // K1
498 //
499
500 RHS(lev,time,
501 density_k1,momentum_k1,energy_k1,
502 density_old, momentum_old, energy_old);
503
504 //
505 // K2
506 //
507
508 // [state] = [old] + (dt/2)[k1]
509 amrex::MultiFab::LinComb(density_st, 1.0, density_old, 0, dt/2.0, density_k1, 0, 0, 1, 0);
510 amrex::MultiFab::LinComb(momentum_st, 1.0, momentum_old, 0, dt/2.0, momentum_k1, 0, 0, 2, 0);
511 amrex::MultiFab::LinComb(energy_st, 1.0, energy_old, 0, dt/2.0, energy_k1, 0, 0, 1, 0);
512
513 density_bc->FillBoundary(density_st, 0, 1, time, 0); density_st.FillBoundary(true);
514 momentum_bc->FillBoundary(momentum_st, 0, 2, time, 0); momentum_st.FillBoundary(true);
515 energy_bc->FillBoundary(energy_st,0,1,time,0); energy_st.FillBoundary(true);
516
517
518 RHS(lev,time,
519 density_k2, momentum_k2, energy_k2,
520 density_st, momentum_st, energy_st);
521
522 //
523 // K3
524 //
525
526 // [state] = [old] + (dt/2)[k2]
527 amrex::MultiFab::LinComb(density_st, 1.0, density_old, 0, dt/2.0, density_k2, 0, 0, 1, 0);
528 amrex::MultiFab::LinComb(momentum_st, 1.0, momentum_old, 0, dt/2.0, momentum_k2, 0, 0, 2, 0);
529 amrex::MultiFab::LinComb(energy_st, 1.0, energy_old, 0, dt/2.0, energy_k2, 0, 0, 1, 0);
530
531 density_bc->FillBoundary(density_st, 0, 1, time, 0); density_st.FillBoundary(true);
532 momentum_bc->FillBoundary(momentum_st, 0, 2, time, 0); momentum_st.FillBoundary(true);
533 energy_bc->FillBoundary(energy_st,0,1,time,0); energy_st.FillBoundary(true);
534
535 RHS(lev,time,
536 density_k3, momentum_k3, energy_k3,
537 density_st, momentum_st, energy_st);
538
539 //
540 // K4
541 //
542
543 // [state] = [old] + (dt/2)[k3]
544 amrex::MultiFab::LinComb(density_st, 1.0, density_old, 0, dt, density_k3, 0, 0, 1, 0);
545 amrex::MultiFab::LinComb(momentum_st, 1.0, momentum_old, 0, dt, momentum_k3, 0, 0, 2, 0);
546 amrex::MultiFab::LinComb(energy_st, 1.0, energy_old, 0, dt, energy_k3, 0, 0, 1, 0);
547
548 density_bc-> FillBoundary(density_st, 0, 1, time, 0); density_st.FillBoundary(true);
549 momentum_bc->FillBoundary(momentum_st, 0, 2, time, 0); momentum_st.FillBoundary(true);
550 energy_bc-> FillBoundary(energy_st, 0, 1, time, 0); energy_st.FillBoundary(true);
551
552
553 RHS(lev,time,
554 density_k4, momentum_k4, energy_k4,
555 density_st, momentum_st, energy_st);
556
557
558 // [new] = [old] + (dt/6)k1
559 amrex::MultiFab::LinComb(density_new, 1.0, density_old, 0, (dt/6.0), density_k1, 0, 0, 1, 0);
560 amrex::MultiFab::LinComb(momentum_new, 1.0, momentum_old, 0, (dt/6.0), momentum_k1, 0, 0, 2, 0);
561 amrex::MultiFab::LinComb(energy_new, 1.0, energy_old, 0, (dt/6.0), energy_k1, 0, 0, 1, 0);
562
563 // [new] += (2 dt/6)k2
564 amrex::MultiFab::Saxpy(density_new, (dt/3.0), density_k2, 0, 0, 1, 0);
565 amrex::MultiFab::Saxpy(momentum_new, (dt/3.0), momentum_k2, 0, 0, 2, 0);
566 amrex::MultiFab::Saxpy(energy_new, (dt/3.0), energy_k2, 0, 0, 1, 0);
567
568 // [new] += (2 dt/6)k3
569 amrex::MultiFab::Saxpy(density_new, (dt/3.0), density_k3, 0, 0, 1, 0);
570 amrex::MultiFab::Saxpy(momentum_new, (dt/3.0), momentum_k3, 0, 0, 2, 0);
571 amrex::MultiFab::Saxpy(energy_new, (dt/3.0), energy_k3, 0, 0, 1, 0);
572
573 // [new] += (dt/6)k4
574 amrex::MultiFab::Saxpy(density_new, (dt/6.0), density_k4, 0, 0, 1, 0);
575 amrex::MultiFab::Saxpy(momentum_new, (dt/6.0), momentum_k4, 0, 0, 2, 0);
576 amrex::MultiFab::Saxpy(energy_new, (dt/6.0), energy_k4, 0, 0, 1, 0);
577
578 }
579
580
581 for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
582 {
583 const amrex::Box& bx = mfi.validbox();
584 const Set::Scalar* DX = geom[lev].CellSize();
585
586 Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
587 Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
588 Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
589 Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
590
591 Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
592 Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
593 Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
594
595 Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
596
597 Set::Patch<Set::Scalar> u = velocity_mf.Patch(lev,mfi);
598 Set::Patch<Set::Scalar> Source = Source_mf.Patch(lev,mfi);
599
600 Set::Scalar *dt_max_handle = &dt_max;
601
602 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
603 {
604 if (eta(i,j,k) < cutoff)
605 {
606 rho_new(i,j,k,0) = rho_solid(i,j,k,0);
607 M_new(i,j,k,0) = M_solid(i,j,k,0);
608 M_new(i,j,k,1) = M_solid(i,j,k,1);
609 E_new(i,j,k,0) = E_solid(i,j,k,0);
610 }
611
612 Set::Matrix gradu = Numeric::Gradient(u, i, j, k, DX);
613 omega(i, j, k) = eta(i, j, k) * (gradu(1,0) - gradu(0,1));
614
615 if (dynamictimestep.on)
616 {
617 *dt_max_handle = std::fabs(cfl * DX[0] / (u(i,j,k,0)*eta(i,j,k) + small));
618 *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl * DX[1] / (u(i,j,k,1)*eta(i,j,k) + small)));
619 *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[0]*DX[0] / (Source(i,j,k,1)+small)));
620 *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[1]*DX[1] / (Source(i,j,k,2)+small)));
621 }
622 });
623 }
624
625
626 if (dynamictimestep.on)
627 {
628 this->DynamicTimestep_SyncTimeStep(lev,dt_max);
629 }
630
631}//end Advance
632
633
634void Hydro::RHS(int lev, Set::Scalar /*time*/,
635 amrex::MultiFab &rho_rhs_mf,
636 amrex::MultiFab &M_rhs_mf,
637 amrex::MultiFab &E_rhs_mf,
638 const amrex::MultiFab &rho_mf,
639 const amrex::MultiFab &M_mf,
640 const amrex::MultiFab &E_mf)
641{
642
643 for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
644 {
645 const amrex::Box& bx = mfi.growntilebox();
646 amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
647
648 Set::Patch<const Set::Scalar> rho = rho_mf.array(mfi); // density
649 Set::Patch<const Set::Scalar> M = M_mf.array(mfi); // momentum
650 Set::Patch<const Set::Scalar> E = E_mf.array(mfi); // total energy (internal energy + kinetic energy) per unit volume (E/rho = e + 0.5*v^2)
651
652 Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
653 Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
654 Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
655
656 Set::Patch<Set::Scalar> v = velocity_mf.Patch(lev,mfi);
657 Set::Patch<Set::Scalar> p = pressure_mf.Patch(lev,mfi);
658
659 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
660 {
661
662 Set::Scalar etarho_fluid = rho(i,j,k) - (1.-eta(i,j,k)) * rho_solid(i,j,k);
663 Set::Scalar etaE_fluid = E(i,j,k) - (1.-eta(i,j,k)) * E_solid(i,j,k);
664
665 Set::Vector etaM_fluid( M(i,j,k,0) - (1.-eta(i,j,k)) * M_solid(i,j,k,0),
666 M(i,j,k,1) - (1.-eta(i,j,k)) * M_solid(i,j,k,1) );
667
668 //THESE ARE FLUID VELOCITY AND PRESSURE
669
670 v(i,j,k,0) = etaM_fluid(0) / (etarho_fluid + small);
671 v(i,j,k,1) = etaM_fluid(1) / (etarho_fluid + small);
672 p(i,j,k) = (etaE_fluid / (eta(i, j, k) + small) - 0.5 * (etaM_fluid(0)*etaM_fluid(0) + etaM_fluid(1)*etaM_fluid(1)) / (etarho_fluid + small)) * ((gamma - 1.0) / (eta(i, j, k) + small))-pref;
673 });
674 }
675
676 const Set::Scalar* DX = geom[lev].CellSize();
677 amrex::Box domain = geom[lev].Domain();
678
679 for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
680 {
681 const amrex::Box& bx = mfi.validbox();
682
683 // Inputs
684 Set::Patch<const Set::Scalar> rho = rho_mf.array(mfi);
685 Set::Patch<const Set::Scalar> E = E_mf.array(mfi);
686 Set::Patch<const Set::Scalar> M = M_mf.array(mfi);
687
688 // Outputs
689 Set::Patch<Set::Scalar> rho_rhs = rho_rhs_mf.array(mfi);
690 Set::Patch<Set::Scalar> M_rhs = M_rhs_mf.array(mfi);
691 Set::Patch<Set::Scalar> E_rhs = E_rhs_mf.array(mfi);
692
693
694 // Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
695 // Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
696 // Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
697
698 Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
699 Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
700 Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
701
702 Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
703
704 Set::Patch<const Set::Scalar> eta = eta_old_mf.Patch(lev,mfi);
705 Set::Patch<const Set::Scalar> etadot = etadot_mf.Patch(lev,mfi);
706 Set::Patch<const Set::Scalar> velocity = velocity_mf.Patch(lev,mfi);
707
708 Set::Patch<const Set::Scalar> m0 = m0_mf.Patch(lev,mfi);
709 Set::Patch<const Set::Scalar> q = q_mf.Patch(lev,mfi);
710 Set::Patch<const Set::Scalar> _u0 = u0_mf.Patch(lev,mfi);
711
712 amrex::Array4<Set::Scalar> const& Source = (*Source_mf[lev]).array(mfi);
713
714 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
715 {
716 auto sten = Numeric::GetStencil(i, j, k, domain);
717
718 //Diffuse Sources
719 Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
720 Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
721 Set::Matrix hess_eta = Numeric::Hessian(eta, i, j, k, 0, DX);
722
723 Set::Vector u = Set::Vector(velocity(i, j, k, 0), velocity(i, j, k, 1));
724 Set::Vector u0 = Set::Vector(_u0(i, j, k, 0), _u0(i, j, k, 1));
725
726 Set::Matrix gradM = Numeric::Gradient(M, i, j, k, DX);
727 Set::Vector gradrho = Numeric::Gradient(rho,i,j,k,0,DX);
728 Set::Matrix hess_rho = Numeric::Hessian(rho,i,j,k,0,DX,sten);
729 Set::Matrix gradu = (gradM - u*gradrho.transpose()) / rho(i,j,k);
730
731 Set::Vector q0 = Set::Vector(q(i,j,k,0),q(i,j,k,1));
732
733 Set::Scalar mdot0 = -m0(i,j,k) * grad_eta_mag;
734 Set::Vector Pdot0 = Set::Vector::Zero();
735 Set::Scalar qdot0 = q0.dot(grad_eta);
736
737 // sten is necessary here because sometimes corner ghost
738 // cells don't get filled
739 Set::Matrix3 hess_M = Numeric::Hessian(M,i,j,k,DX);
740 Set::Matrix3 hess_u = Set::Matrix3::Zero();
741 for (int p = 0; p < 2; p++)
742 for (int q = 0; q < 2; q++)
743 for (int r = 0; r < 2; r++)
744 {
745 hess_u(r,p,q) =
746 (hess_M(r,p,q) - gradu(r,q)*gradrho(p) - gradu(r,p)*gradrho(q) - u(r)*hess_rho(p,q))
747 / rho(i,j,k);
748 }
749
750 Set::Vector Ldot0 = Set::Vector::Zero();
751 Set::Vector div_tau = Set::Vector::Zero();
752 for (int p = 0; p<2; p++)
753 for (int q = 0; q<2; q++)
754 for (int r = 0; r<2; r++)
755 for (int s = 0; s<2; s++)
756 {
757 Set::Scalar Mpqrs = 0.0;
758 if (p==r && q==s) Mpqrs += 0.5 * mu;
759
760 Ldot0(p) += 0.5*Mpqrs * (u(r) - u0(r)) * hess_eta(q, s);
761 div_tau(p) += 2.0*Mpqrs * hess_u(r,s,q);
762 }
763
764 Source(i,j, k, 0) = mdot0;
765 Source(i,j, k, 1) = Pdot0(0) - Ldot0(0);
766 Source(i,j, k, 2) = Pdot0(1) - Ldot0(1);
767 Source(i,j, k, 3) = qdot0;// - Ldot0(0)*v(i,j,k,0) - Ldot0(1)*v(i,j,k,1);
768
769 // Lagrange terms to enforce no-penetration
770 Source(i,j,k,1) -= lagrange*(u-u0).dot(grad_eta)*grad_eta(0);
771 Source(i,j,k,2) -= lagrange*(u-u0).dot(grad_eta)*grad_eta(1);
772
773 //Godunov flux
774 //states of total fields
775 const int X = 0, Y = 1;
776 Solver::Local::Riemann::State state_xlo(rho, M, E, i-1, j, k, X);
777 Solver::Local::Riemann::State state_x (rho, M, E, i , j, k, X);
778 Solver::Local::Riemann::State state_xhi(rho, M, E, i+1, j, k, X);
779
780 Solver::Local::Riemann::State state_ylo(rho, M, E, i, j-1, k, Y);
781 Solver::Local::Riemann::State state_y (rho, M, E, i, j , k, Y);
782 Solver::Local::Riemann::State state_yhi(rho, M, E, i, j+1, k, Y);
783
784 //states of solid fields
785 Solver::Local::Riemann::State state_xlo_solid(rho_solid, M_solid, E_solid, i-1, j, k, X);
786 Solver::Local::Riemann::State state_x_solid (rho_solid, M_solid, E_solid, i , j, k, X);
787 Solver::Local::Riemann::State state_xhi_solid(rho_solid, M_solid, E_solid, i+1, j, k, X);
788
789 Solver::Local::Riemann::State state_ylo_solid(rho_solid, M_solid, E_solid, i, j-1, k, Y);
790 Solver::Local::Riemann::State state_y_solid (rho_solid, M_solid, E_solid, i, j , k, Y);
791 Solver::Local::Riemann::State state_yhi_solid(rho_solid, M_solid, E_solid, i, j+1, k, Y);
792
793 Solver::Local::Riemann::State state_xlo_fluid = (state_xlo - (1.0 - eta(i-1,j,k))*state_xlo_solid) / (eta(i-1,j,k) + small);
794 Solver::Local::Riemann::State state_x_fluid = (state_x - (1.0 - eta(i,j,k) )*state_x_solid ) / (eta(i,j,k) + small);
795 Solver::Local::Riemann::State state_xhi_fluid = (state_xhi - (1.0 - eta(i+1,j,k))*state_xhi_solid) / (eta(i+1,j,k) + small);
796
797 Solver::Local::Riemann::State state_ylo_fluid = (state_ylo - (1.0 - eta(i,j-1,k))*state_ylo_solid) / (eta(i,j-1,k) + small);
798 Solver::Local::Riemann::State state_y_fluid = (state_y - (1.0 - eta(i,j,k) )*state_y_solid ) / (eta(i,j,k) + small);
799 Solver::Local::Riemann::State state_yhi_fluid = (state_yhi - (1.0 - eta(i,j+1,k))*state_yhi_solid) / (eta(i,j+1,k) + small);
800
801 Solver::Local::Riemann::Flux flux_xlo, flux_ylo, flux_xhi, flux_yhi;
802
803 try
804 {
805 //lo interface fluxes
806 flux_xlo = riemannsolver->Solve(state_xlo_fluid, state_x_fluid, gamma, pref, small) * eta(i,j,k);
807 flux_ylo = riemannsolver->Solve(state_ylo_fluid, state_y_fluid, gamma, pref, small) * eta(i,j,k);
808
809 //hi interface fluxes
810 flux_xhi = riemannsolver->Solve(state_x_fluid, state_xhi_fluid, gamma, pref, small) * eta(i,j,k);
811 flux_yhi = riemannsolver->Solve(state_y_fluid, state_yhi_fluid, gamma, pref, small) * eta(i,j,k);
812 }
813 catch(...)
814 {
815 Util::ParallelMessage(INFO,"lev=",lev);
816 Util::ParallelMessage(INFO,"i=",i,"j=",j);
817 Util::Abort(INFO);
818 }
819
820
821 Set::Scalar drhof_dt =
822 (flux_xlo.mass - flux_xhi.mass) / DX[0] +
823 (flux_ylo.mass - flux_yhi.mass) / DX[1] +
824 Source(i, j, k, 0);
825
826 rho_rhs(i,j,k) =
827 // rho_new(i, j, k) = rho(i, j, k) +
828 //(
829 drhof_dt +
830 // todo add drhos_dt term if want time-evolving rhos
831 etadot(i,j,k) * (rho(i,j,k) - rho_solid(i,j,k)) / (eta(i,j,k) + small)
832 // ) * dt;
833 ;
834
835
836
837 Set::Scalar dMxf_dt =
838 (flux_xlo.momentum_normal - flux_xhi.momentum_normal ) / DX[0] +
839 (flux_ylo.momentum_tangent - flux_yhi.momentum_tangent) / DX[1] +
840 div_tau(0) * eta(i,j,k) +
841 Source(i, j, k, 1);
842
843 M_rhs(i,j,k,0) =
844 //M_new(i, j, k, 0) = M(i, j, k, 0) +
845 // (
846 dMxf_dt +
847 // todo add dMs_dt term if want time-evolving Ms
848 etadot(i,j,k)*(M(i,j,k,0) - M_solid(i,j,k,0)) / (eta(i,j,k) + small)
849 // ) * dt;
850 ;
851
852 Set::Scalar dMyf_dt =
853 (flux_xlo.momentum_tangent - flux_xhi.momentum_tangent) / DX[0] +
854 (flux_ylo.momentum_normal - flux_yhi.momentum_normal ) / DX[1] +
855 div_tau(1) * eta(i,j,k) +
856 Source(i, j, k, 2);
857
858 M_rhs(i,j,k,1) =
859 //M_new(i, j, k, 1) = M(i, j, k, 1) +
860 //(
861 dMyf_dt +
862 // todo add dMs_dt term if want time-evolving Ms
863 etadot(i,j,k)*(M(i,j,k,1) - M_solid(i,j,k,1)) / (eta(i,j,k)+small)
864 // )*dt;
865 ;
866
867 Set::Scalar dEf_dt =
868 (flux_xlo.energy - flux_xhi.energy) / DX[0] +
869 (flux_ylo.energy - flux_yhi.energy) / DX[1] +
870 Source(i, j, k, 3);
871
872 E_rhs(i,j,k) =
873 // E_new(i, j, k) = E(i, j, k) +
874 // (
875 dEf_dt +
876 // todo add dEs_dt term if want time-evolving Es
877 etadot(i,j,k)*(E(i,j,k) - E_solid(i,j,k)) / (eta(i,j,k)+small)
878 // ) * dt;
879 ;
880
881#ifdef AMREX_DEBUG
882 if ((rho_rhs(i,j,k) != rho_rhs(i,j,k)) ||
883 (M_rhs(i,j,k,0) != M_rhs(i,j,k,0)) ||
884 (M_rhs(i,j,k,1) != M_rhs(i,j,k,1)) ||
885 (E_rhs(i,j,k) != E_rhs(i,j,k)))
886 {
887 Util::ParallelMessage(INFO,"rho_rhs=",rho_rhs(i,j,k));
888 Util::ParallelMessage(INFO,"Mx_rhs=",M_rhs(i,j,k,0));
889 Util::ParallelMessage(INFO,"Mx_rhs=",M_rhs(i,j,k,1));
890 Util::ParallelMessage(INFO,"E_rhs=",E_rhs(i,j,k));
891
892 Util::ParallelMessage(INFO,"lev=",lev);
893 Util::ParallelMessage(INFO,"i=",i," j=",j);
894 Util::ParallelMessage(INFO,"drhof_dt ",drhof_dt); // dies
895 Util::ParallelMessage(INFO,"flux_xlo.mass ",flux_xlo.mass);
896 Util::ParallelMessage(INFO,"flux_xhi.mass ",flux_xhi.mass); // dies, depends on state_xx, state_xhi, state_x_solid, state_xhi_solid, gamma, eta, pref, small
897 Util::ParallelMessage(INFO,"flux_ylo.mass ",flux_ylo.mass);
898 Util::ParallelMessage(INFO,"flux_xhi.mass ",flux_yhi.mass);
899 Util::ParallelMessage(INFO,"eta ",eta(i,j,k));
900 Util::ParallelMessage(INFO,"etadot ",etadot(i,j,k));
901 Util::ParallelMessage(INFO,"Source ",Source(i,j,k,0));
902 Util::ParallelMessage(INFO,"state_x ",state_x); // <<<<
903 Util::ParallelMessage(INFO,"state_y ",state_y);
904 Util::ParallelMessage(INFO,"state_x_solid ",state_x_solid); // <<<<
905 Util::ParallelMessage(INFO,"state_y_solid ",state_y_solid);
906 Util::ParallelMessage(INFO,"state_xhi ",state_xhi); // <<<<
907 Util::ParallelMessage(INFO,"state_yhi ",state_yhi);
908 Util::ParallelMessage(INFO,"state_xhi_solid ",state_xhi_solid);
909 Util::ParallelMessage(INFO,"state_yhi_solids ",state_yhi_solid);
910 Util::ParallelMessage(INFO,"state_xlo ",state_xlo);
911 Util::ParallelMessage(INFO,"state_ylo ",state_ylo);
912 Util::ParallelMessage(INFO,"state_xlo_solid ",state_xlo_solid);
913 Util::ParallelMessage(INFO,"state_ylo_solid ",state_ylo_solid);
914
915 Util::ParallelMessage(INFO,"Mx_solid ",M_solid(i,j,k,0));
916 Util::ParallelMessage(INFO,"My_solid ",M_solid(i,j,k,1));
917 Util::ParallelMessage(INFO,"small ",small);
918 Util::ParallelMessage(INFO,"Mx ",M(i,j,k,0));
919 Util::ParallelMessage(INFO,"My ",M(i,j,k,1));
920 Util::ParallelMessage(INFO,"dMx/dt ",dMxf_dt);
921 Util::ParallelMessage(INFO,"dMy/dt ",dMyf_dt);
922
923
924 Util::Message(INFO,flux_xlo.momentum_tangent);
925 Util::Message(INFO,flux_xhi.momentum_tangent);
926 Util::Message(INFO,DX[0]);
927 Util::Message(INFO,flux_ylo.momentum_normal);
928 Util::Message(INFO,flux_yhi.momentum_normal);
929 Util::Message(INFO,DX[1]);
930 Util::Message(INFO,div_tau);
931 Util::Message(INFO,Source(i, j, k, 2));
932
933 Util::Message(INFO,hess_eta);
934 Util::Message(INFO,velocity(i,j,k,0));
935 Util::Message(INFO,velocity(i,j,k,1));
936
937 Util::Exception(INFO);
938 }
939#endif
940
941
942
943 // todo - may need to move this for higher order schemes...
944 omega(i, j, k) = eta(i, j, k) * (gradu(1,0) - gradu(0,1));
945 });
946 }
947}
948
949void Hydro::Regrid(int lev, Set::Scalar /* time */)
950{
951 BL_PROFILE("Integrator::Hydro::Regrid");
952 Source_mf[lev]->setVal(0.0);
953 if (lev < finest_level) return;
954
955 Util::Message(INFO, "Regridding on level", lev);
956}//end regrid
957
958//void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, Set::Scalar time, int ngrow)
959void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar, int)
960{
961 BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
962
963 const Set::Scalar* DX = geom[lev].CellSize();
964 Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
965
966 // Eta criterion for refinement
967 for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi) {
968 const amrex::Box& bx = mfi.tilebox();
969 amrex::Array4<char> const& tags = a_tags.array(mfi);
970 amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
971
972 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
973 Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
974 if (grad_eta.lpNorm<2>() * dr * 2 > eta_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
975 });
976 }
977
978 // Vorticity criterion for refinement
979 for (amrex::MFIter mfi(*vorticity_mf[lev], true); mfi.isValid(); ++mfi) {
980 const amrex::Box& bx = mfi.tilebox();
981 amrex::Array4<char> const& tags = a_tags.array(mfi);
982 amrex::Array4<const Set::Scalar> const& omega = (*vorticity_mf[lev]).array(mfi);
983
984 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
985 auto sten = Numeric::GetStencil(i, j, k, bx);
986 Set::Vector grad_omega = Numeric::Gradient(omega, i, j, k, 0, DX, sten);
987 if (grad_omega.lpNorm<2>() * dr * 2 > omega_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
988 });
989 }
990
991 // Gradu criterion for refinement
992 for (amrex::MFIter mfi(*velocity_mf[lev], true); mfi.isValid(); ++mfi) {
993 const amrex::Box& bx = mfi.tilebox();
994 amrex::Array4<char> const& tags = a_tags.array(mfi);
995 amrex::Array4<const Set::Scalar> const& v = (*velocity_mf[lev]).array(mfi);
996
997 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
998 auto sten = Numeric::GetStencil(i, j, k, bx);
999 Set::Matrix grad_u = Numeric::Gradient(v, i, j, k, DX, sten);
1000 if (grad_u.lpNorm<2>() * dr * 2 > gradu_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
1001 });
1002 }
1003
1004 // Pressure criterion for refinement
1005 for (amrex::MFIter mfi(*pressure_mf[lev], true); mfi.isValid(); ++mfi) {
1006 const amrex::Box& bx = mfi.tilebox();
1007 amrex::Array4<char> const& tags = a_tags.array(mfi);
1008 amrex::Array4<const Set::Scalar> const& p = (*pressure_mf[lev]).array(mfi);
1009
1010 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
1011 auto sten = Numeric::GetStencil(i, j, k, bx);
1012 Set::Vector grad_p = Numeric::Gradient(p, i, j, k, 0, DX, sten);
1013 if (grad_p.lpNorm<2>() * dr * 2 > p_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
1014 });
1015 }
1016
1017 // Density criterion for refinement
1018 for (amrex::MFIter mfi(*density_mf[lev], true); mfi.isValid(); ++mfi) {
1019 const amrex::Box& bx = mfi.tilebox();
1020 amrex::Array4<char> const& tags = a_tags.array(mfi);
1021 amrex::Array4<const Set::Scalar> const& rho = (*density_mf[lev]).array(mfi);
1022
1023 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
1024 auto sten = Numeric::GetStencil(i, j, k, bx);
1025 Set::Vector grad_rho = Numeric::Gradient(rho, i, j, k, 0, DX, sten);
1026 if (grad_rho.lpNorm<2>() * dr * 2 > rho_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
1027 });
1028 }
1029
1030}
1031
1032}
1033
1034
1035#endif
X
#define X(name)
pp_query_required
#define pp_query_required(...)
Definition ParmParse.H:101
pp_query_default
#define pp_query_default(...)
Definition ParmParse.H:102
pp_forbid
#define pp_forbid(...)
Definition ParmParse.H:110
pp_queryclass
#define pp_queryclass(...)
Definition ParmParse.H:109
TEST
#define TEST(x)
Definition Util.H:22
INFO
#define INFO
Definition Util.H:21
IC::BMP
Definition BMP.H:22
IC::Laminate
Initialize Laminates in a matrix.
Definition Laminate.H:16
IC::PNG
Definition PNG.H:26
IO::ParmParse::forbid
void forbid(std::string name, std::string explanation)
Definition ParmParse.H:160
IO::ParmParse::AnyUnusedInputs
int AnyUnusedInputs(bool inscopeonly=true, bool verbose=false)
Definition ParmParse.H:869
IO::ParmParse::select_default
void select_default(std::string name, PTRTYPE *&ic_eta, Args &&... args)
Definition ParmParse.H:1069
IO::ParmParse::AllUnusedInputs
static int AllUnusedInputs()
Definition ParmParse.H:907
IO::ParmParse::query_default
int query_default(std::string name, T &value, T defaultvalue)
Definition ParmParse.H:293
IO::ParmParse::query_validate
int query_validate(std::string name, int &value, std::vector< int > possibleintvals)
Definition ParmParse.H:336
Set::Matrix3::Zero
static Matrix3 Zero()
Definition Matrix3.H:42
Solver::Local::Riemann::Roe
Roe Riemann Solver based on Gas Dynamics - Culbert B. Laney.
Definition Roe.H:30
Integrator
Collection of numerical integrator objects.
Definition AllenCahn.H:42
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::Hessian
AMREX_FORCE_INLINE Set::Matrix Hessian(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:966
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
Set::Scalar
amrex::Real Scalar
Definition Base.H:18
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:19
Set::Matrix
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, AMREX_SPACEDIM > Matrix
Definition Base.H:22
Util::ParallelMessage
void ParallelMessage(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:182
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:225
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:55
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:199
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:126
Util::Exception
void Exception(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:223
Solver::Local::Riemann::Flux::momentum_normal
Set::Scalar momentum_normal
Definition Riemann.H:92
Solver::Local::Riemann::Flux::momentum_tangent
Set::Scalar momentum_tangent
Definition Riemann.H:93
Generated by doxygen 1.9.8