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