Line data Source code
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 :
17 : namespace Integrator
18 : {
19 :
20 6 : Hydro::Hydro(IO::ParmParse& pp) : Hydro()
21 : {
22 18 : pp_queryclass(*this);
23 6 : }
24 :
25 : void
26 6 : Hydro::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 6 : std::string scheme_str;
37 : // time integration scheme to use
38 42 : pp.query_validate("scheme",scheme_str, {"forwardeuler","ssprk3","rk4"});
39 6 : if (scheme_str == "forwardeuler") value.scheme = IntegrationScheme::ForwardEuler;
40 2 : else if (scheme_str == "ssprk3") value.scheme = IntegrationScheme::SSPRK3;
41 1 : else if (scheme_str == "rk4") value.scheme = IntegrationScheme::RK4;
42 :
43 :
44 : // eta-based refinement
45 36 : pp.query_default("eta_refinement_criterion", value.eta_refinement_criterion , 0.01);
46 : // vorticity-based refinement
47 36 : pp.query_default("omega_refinement_criterion", value.omega_refinement_criterion, 0.01);
48 : // velocity gradient-based refinement
49 36 : pp.query_default("gradu_refinement_criterion", value.gradu_refinement_criterion, 0.01);
50 : // pressure-based refinement
51 36 : pp.query_default("p_refinement_criterion", value.p_refinement_criterion, 1e100);
52 : // density-based refinement
53 36 : pp.query_default("rho_refinement_criterion", value.rho_refinement_criterion, 1e100);
54 :
55 36 : pp_query_required("gamma", value.gamma); // gamma for gamma law
56 36 : pp_query_required("cfl", value.cfl); // cfl condition
57 36 : pp_query_default("cfl_v", value.cfl_v,1E100); // cfl condition
58 36 : pp_query_required("mu", value.mu); // linear viscosity coefficient
59 48 : pp_forbid("Lfactor","replaced with mu");
60 : //pp_query_default("Lfactor", value.Lfactor,1.0); // (to be removed) test factor for viscous source
61 48 : pp_forbid("Pfactor","replaced with mu");
62 : //pp_query_default("Pfactor", value.Pfactor,1.0); // (to be removed) test factor for viscous source
63 36 : pp_query_default("pref", value.pref,1.0); // reference pressure for Roe solver
64 :
65 48 : pp_forbid("rho.bc","--> density.bc");
66 48 : pp_forbid("p.bc","--> pressure.bc");
67 48 : pp_forbid("v.bc", "--> velocity.bc");
68 48 : pp_forbid("pressure.bc","--> energy.bc");
69 42 : pp_forbid("velocity.bc","--> momentum.bc");
70 :
71 : // Boundary condition for density
72 12 : pp.select_default<BC::Constant,BC::Expression>("density.bc",value.density_bc,1);
73 : // Boundary condition for energy
74 12 : pp.select_default<BC::Constant,BC::Expression>("energy.bc",value.energy_bc,1);
75 : // Boundary condition for momentum
76 12 : pp.select_default<BC::Constant,BC::Expression>("momentum.bc",value.momentum_bc,2);
77 : // Boundary condition for phase field order parameter
78 18 : pp.select_default<BC::Constant,BC::Expression>("pf.eta.bc",value.eta_bc,1);
79 :
80 36 : pp_query_default("small",value.small,1E-8); // small regularization value
81 36 : pp_query_default("cutoff",value.cutoff,-1E100); // cutoff value
82 36 : pp_query_default("lagrange",value.lagrange,0.0); // lagrange no-penetration factor
83 :
84 42 : pp_forbid("roefix","--> solver.roe.entropy_fix"); // Roe solver entropy fix
85 :
86 6 : }
87 : // Register FabFields:
88 : {
89 6 : int nghost = 1;
90 :
91 18 : value.RegisterNewFab(value.eta_mf, value.eta_bc, 1, nghost, "eta", true, true);
92 18 : value.RegisterNewFab(value.eta_old_mf, value.eta_bc, 1, nghost, "eta_old", true, true);
93 18 : value.RegisterNewFab(value.etadot_mf, value.eta_bc, 1, nghost, "etadot", true, false);
94 :
95 18 : value.RegisterNewFab(value.density_mf, value.density_bc, 1, nghost, "density", true , true);
96 18 : value.RegisterNewFab(value.density_old_mf, value.density_bc, 1, nghost, "density_old", false, true);
97 :
98 18 : value.RegisterNewFab(value.energy_mf, value.energy_bc, 1, nghost, "energy", true ,true);
99 18 : value.RegisterNewFab(value.energy_old_mf, value.energy_bc, 1, nghost, "energy_old" , false, true);
100 :
101 24 : value.RegisterNewFab(value.momentum_mf, value.momentum_bc, 2, nghost, "momentum", true ,true, {"x","y"});
102 18 : value.RegisterNewFab(value.momentum_old_mf, value.momentum_bc, 2, nghost, "momentum_old", false, true);
103 :
104 18 : value.RegisterNewFab(value.pressure_mf, &value.bc_nothing, 1, nghost, "pressure", true, false);
105 24 : value.RegisterNewFab(value.velocity_mf, &value.bc_nothing, 2, nghost, "velocity", true, false,{"x","y"});
106 18 : value.RegisterNewFab(value.vorticity_mf, &value.bc_nothing, 1, nghost, "vorticity", true, false);
107 :
108 18 : value.RegisterNewFab(value.m0_mf, &value.bc_nothing, 1, 0, "m0", true, false);
109 24 : value.RegisterNewFab(value.u0_mf, &value.bc_nothing, 2, 0, "u0", true, false, {"x","y"});
110 24 : value.RegisterNewFab(value.q_mf, &value.bc_nothing, 2, 0, "q", true, false, {"x","y"});
111 :
112 24 : value.RegisterNewFab(value.solid.momentum_mf, &value.neumann_bc_D, 2, nghost, "solid.momentum", true, false, {"x","y"});
113 18 : value.RegisterNewFab(value.solid.density_mf, &value.neumann_bc_1, 1, nghost, "solid.density", true, false);
114 18 : value.RegisterNewFab(value.solid.energy_mf, &value.neumann_bc_1, 1, nghost, "solid.energy", true, false);
115 :
116 18 : value.RegisterNewFab(value.Source_mf, &value.bc_nothing, 4, 0, "Source", true, false);
117 : }
118 :
119 48 : pp_forbid("Velocity.ic.type", "--> velocity.ic.type");
120 48 : pp_forbid("Pressure.ic", "--> pressure.ic");
121 48 : pp_forbid("SolidMomentum.ic", "--> solid.momentum.ic");
122 48 : pp_forbid("SolidDensity.ic.type", "--> solid.density.ic.type");
123 48 : pp_forbid("SolidEnergy.ic.type", "--> solid.energy.ic.type");
124 48 : pp_forbid("Density.ic.type", "--> density.ic.type");
125 48 : pp_forbid("rho_injected.ic.type","no longer using rho_injected use m0 instead");
126 42 : pp.forbid("mdot.ic.type", "replace mdot with u0");
127 :
128 :
129 : // ORDER PARAMETER
130 :
131 : // eta initial condition
132 12 : 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 12 : pp.select_default<IC::Constant,IC::Expression>("velocity.ic",value.velocity_ic,value.geom);
138 : // solid pressure initial condition
139 12 : pp.select_default<IC::Constant,IC::Expression>("pressure.ic",value.pressure_ic,value.geom);
140 : // density initial condition type
141 12 : 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 12 : pp.select_default<IC::Constant,IC::Expression>("solid.momentum.ic",value.solid.momentum_ic,value.geom);
148 : // solid density initial condition
149 12 : pp.select_default<IC::Constant,IC::Expression>("solid.density.ic",value.solid.density_ic,value.geom);
150 : // solid energy initial condition
151 12 : 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 12 : pp.select_default<IC::Constant,IC::Expression>("m0.ic",value.ic_m0,value.geom);
158 : // diffuse boundary prescribed velocity
159 12 : pp.select_default<IC::Constant,IC::Expression>("u0.ic",value.ic_u0,value.geom);
160 : // diffuse boundary prescribed heat flux
161 12 : 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 18 : 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 30 : pp.query_default("allow_unused",allow_unused,false);
174 6 : if (!allow_unused && pp.AnyUnusedInputs(true, false))
175 : {
176 0 : Util::Warning(INFO,"The following inputs were specified but not used:");
177 0 : pp.AllUnusedInputs();
178 0 : Util::Exception(INFO,"Aborting. Specify 'allow_unused=True` to ignore this error.");
179 : }
180 6 : }
181 :
182 :
183 6 : void Hydro::Initialize(int lev)
184 : {
185 : BL_PROFILE("Integrator::Hydro::Initialize");
186 :
187 6 : eta_ic ->Initialize(lev, eta_mf, 0.0);
188 6 : eta_ic ->Initialize(lev, eta_old_mf, 0.0);
189 6 : etadot_mf[lev] ->setVal(0.0);
190 :
191 : //flux_mf[lev] ->setVal(0.0);
192 :
193 6 : velocity_ic ->Initialize(lev, velocity_mf, 0.0);
194 6 : pressure_ic ->Initialize(lev, pressure_mf, 0.0);
195 6 : density_ic ->Initialize(lev, density_mf, 0.0);
196 :
197 6 : density_ic ->Initialize(lev, density_old_mf, 0.0);
198 :
199 :
200 6 : solid.density_ic ->Initialize(lev, solid.density_mf, 0.0);
201 6 : solid.momentum_ic->Initialize(lev, solid.momentum_mf, 0.0);
202 6 : solid.energy_ic ->Initialize(lev, solid.energy_mf, 0.0);
203 :
204 6 : ic_m0 ->Initialize(lev, m0_mf, 0.0);
205 6 : ic_u0 ->Initialize(lev, u0_mf, 0.0);
206 6 : ic_q ->Initialize(lev, q_mf, 0.0);
207 :
208 6 : Source_mf[lev] ->setVal(0.0);
209 :
210 6 : Mix(lev);
211 6 : }
212 :
213 6 : void Hydro::Mix(int lev)
214 : {
215 12 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
216 : {
217 6 : const amrex::Box& bx = mfi.growntilebox();
218 :
219 6 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
220 :
221 6 : Set::Patch<const Set::Scalar> v = velocity_mf.Patch(lev,mfi);
222 6 : Set::Patch<const Set::Scalar> p = pressure_mf.Patch(lev,mfi);
223 6 : Set::Patch<Set::Scalar> rho = density_mf.Patch(lev,mfi);
224 6 : Set::Patch<Set::Scalar> rho_old = density_old_mf.Patch(lev,mfi);
225 6 : Set::Patch<Set::Scalar> M = momentum_mf.Patch(lev,mfi);
226 6 : Set::Patch<Set::Scalar> M_old = momentum_old_mf.Patch(lev,mfi);
227 6 : Set::Patch<Set::Scalar> E = energy_mf.Patch(lev,mfi);
228 6 : Set::Patch<Set::Scalar> E_old = energy_old_mf.Patch(lev,mfi);
229 6 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
230 6 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
231 6 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
232 :
233 :
234 6 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
235 : {
236 80400 : rho(i, j, k) = eta(i, j, k) * rho(i, j, k) + (1.0 - eta(i, j, k)) * rho_solid(i, j, k);
237 32160 : rho_old(i, j, k) = rho(i, j, k);
238 :
239 96480 : 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 96480 : 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 32160 : M_old(i, j, k, 0) = M(i, j, k, 0);
242 32160 : M_old(i, j, k, 1) = M(i, j, k, 1);
243 :
244 32160 : E(i, j, k) =
245 112560 : (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 16080 : +
247 48240 : E_solid(i, j, k) * (1.0 - eta(i, j, k));
248 32160 : E_old(i, j, k) = E(i, j, k);
249 16080 : });
250 6 : }
251 6 : c_max = 0.0;
252 6 : vx_max = 0.0;
253 6 : vy_max = 0.0;
254 6 : }
255 :
256 3650 : void Hydro::UpdateEta(int lev, Set::Scalar time)
257 : {
258 3650 : eta_ic->Initialize(lev, eta_mf, time);
259 3650 : }
260 :
261 3650 : void Hydro::TimeStepBegin(Set::Scalar, int /*iter*/)
262 : {
263 :
264 3650 : }
265 :
266 3650 : void Hydro::TimeStepComplete(Set::Scalar, int lev)
267 : {
268 3650 : if (dynamictimestep.on)
269 0 : Integrator::DynamicTimestep_Update();
270 3650 : 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 :
285 3650 : void Hydro::Advance(int lev, Set::Scalar time, Set::Scalar dt)
286 : {
287 :
288 3650 : std::swap(eta_old_mf, eta_mf);
289 3650 : std::swap(density_old_mf[lev], density_mf[lev]);
290 3650 : std::swap(momentum_old_mf[lev], momentum_mf[lev]);
291 3650 : std::swap(energy_old_mf[lev], energy_mf[lev]);
292 3650 : Set::Scalar dt_max = std::numeric_limits<Set::Scalar>::max();
293 :
294 3650 : UpdateEta(lev, time);
295 :
296 7300 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
297 : {
298 3650 : const amrex::Box& bx = mfi.growntilebox();
299 3650 : amrex::Array4<const Set::Scalar> const& eta_new = (*eta_mf[lev]).array(mfi);
300 3650 : amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
301 3650 : amrex::Array4<Set::Scalar> const& etadot = (*etadot_mf[lev]).array(mfi);
302 3650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
303 : {
304 31951800 : etadot(i, j, k) = (eta_new(i, j, k) - eta(i, j, k)) / dt;
305 10650600 : });
306 3650 : }
307 :
308 :
309 3650 : if (scheme == IntegrationScheme::ForwardEuler) // forward euler
310 : {
311 :
312 3250 : RHS(lev, time,
313 3250 : *density_mf[lev], *momentum_mf[lev], *energy_mf[lev],
314 3250 : *density_old_mf[lev], *momentum_old_mf[lev], *energy_old_mf[lev]);
315 :
316 6500 : for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
317 : {
318 3250 : const amrex::Box& bx = mfi.validbox();
319 :
320 3250 : Set::Patch<const Set::Scalar> rho_rhs = density_mf.Patch(lev,mfi);
321 3250 : Set::Patch<const Set::Scalar> E_rhs = energy_mf.Patch(lev,mfi);
322 3250 : Set::Patch<const Set::Scalar> M_rhs = momentum_mf.Patch(lev,mfi);
323 :
324 3250 : Set::Patch<const Set::Scalar> rho_old = density_old_mf.Patch(lev,mfi);
325 3250 : Set::Patch<const Set::Scalar> E_old = energy_old_mf.Patch(lev,mfi);
326 3250 : Set::Patch<const Set::Scalar> M_old = momentum_old_mf.Patch(lev,mfi);
327 :
328 3250 : Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
329 3250 : Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
330 3250 : Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
331 :
332 :
333 3250 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
334 : {
335 18816000 : rho_new(i, j, k) = rho_old(i, j, k) + dt * rho_rhs(i,j,k);
336 18816000 : M_new(i,j,k,0) = M_old(i,j,k,0) + dt * M_rhs(i,j,k,0);
337 18816000 : M_new(i,j,k,1) = M_old(i,j,k,1) + dt * M_rhs(i,j,k,1);
338 18816000 : E_new(i,j,k) = E_old(i,j,k) + dt * E_rhs(i,j,k);
339 6272000 : });
340 3250 : }
341 : }
342 :
343 :
344 400 : 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 200 : /* */ c2 = 1.0 , a21 = 1.0,
356 200 : /* */ c3 = 0.5, a31 = 0.25, a32 = 0.25,
357 : /* --------------------------------------------- */
358 200 : /* */ b1 = 1./6, b2 = 1./6., b3 = 2./3.;
359 :
360 :
361 200 : const amrex::BoxArray &ba = density_mf[lev]->boxArray();
362 200 : const amrex::DistributionMapping &dm = density_mf[lev]->DistributionMap();
363 200 : const int ng = density_mf[lev]->nGrow();
364 :
365 : // handles to old solution
366 200 : const amrex::MultiFab &density_old = *density_old_mf[lev];
367 200 : const amrex::MultiFab &momentum_old = *momentum_old_mf[lev];
368 200 : const amrex::MultiFab &energy_old = *energy_old_mf[lev];
369 :
370 :
371 : // temporary storage
372 200 : amrex::MultiFab density_k1(ba,dm,1,0), momentum_k1(ba,dm,2,0), energy_k1(ba,dm,1,0);
373 200 : amrex::MultiFab density_k2(ba,dm,1,0), momentum_k2(ba,dm,2,0), energy_k2(ba,dm,1,0);
374 200 : 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 200 : 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 200 : density_temp.ParallelCopyToGhost(*density_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
381 200 : momentum_temp.ParallelCopyToGhost(*momentum_old_mf[lev],0,0,2,amrex::IntVect(1),amrex::IntVect(1));
382 200 : energy_temp.ParallelCopyToGhost(*energy_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
383 :
384 : // handles to new solution
385 200 : amrex::MultiFab &density_new = *density_mf[lev];
386 200 : amrex::MultiFab &momentum_new = *momentum_mf[lev];
387 200 : amrex::MultiFab &energy_new = *energy_mf[lev];
388 :
389 :
390 : //
391 : // Calculate K1
392 : //
393 : // k1 = RHS(t, yold)
394 :
395 200 : RHS(lev,time,
396 : density_k1,momentum_k1,energy_k1,
397 200 : *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 200 : amrex::MultiFab::LinComb(density_temp, 1.0, density_old, 0, dt*a21, density_k1, 0, 0, 1, 0);
405 200 : amrex::MultiFab::LinComb(momentum_temp, 1.0, momentum_old, 0, dt*a21, momentum_k1, 0, 0, 2, 0);
406 200 : amrex::MultiFab::LinComb(energy_temp, 1.0, energy_old, 0, dt*a21, energy_k1, 0, 0, 1, 0);
407 : // fill boundary
408 200 : density_bc ->FillBoundary(density_temp, 0, 1, time, 0); density_temp.FillBoundary(true);
409 200 : momentum_bc->FillBoundary(momentum_temp, 0, 2, time, 0); momentum_temp.FillBoundary(true);
410 200 : energy_bc ->FillBoundary(energy_temp, 0, 1, time, 0); energy_temp.FillBoundary(true);
411 : // k2 = RHS(t + c2*dt, ytemp)
412 200 : 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 200 : amrex::MultiFab::LinComb(density_temp, 1.0, density_old, 0, dt*a31, density_k1, 0, 0, 1, 0);
423 200 : amrex::MultiFab::LinComb(momentum_temp, 1.0, momentum_old, 0, dt*a31, momentum_k1, 0, 0, 2, 0);
424 200 : 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 200 : amrex::MultiFab::Saxpy(density_temp, dt*a32, density_k2, 0, 0, 1, 0);
427 200 : amrex::MultiFab::Saxpy(momentum_temp, dt*a32, momentum_k2, 0, 0, 2, 0);
428 200 : amrex::MultiFab::Saxpy(energy_temp, dt*a32, energy_k2, 0, 0, 1, 0);
429 : // 3. fill boundary
430 200 : density_bc ->FillBoundary(density_temp, 0, 1, time+c2*dt, 0); density_temp.FillBoundary(true);
431 200 : momentum_bc->FillBoundary(momentum_temp, 0, 2, time+c2*dt, 0); momentum_temp.FillBoundary(true);
432 200 : energy_bc ->FillBoundary(energy_temp, 0, 1, time+c2*dt, 0); energy_temp.FillBoundary(true);
433 : // 4. k3 = RHS(t + c3*dt, ytemp)
434 200 : 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 200 : amrex::MultiFab::LinComb(density_new, 1.0, density_old, 0, dt*b1, density_k1, 0, 0, 1, 0);
445 200 : amrex::MultiFab::LinComb(momentum_new, 1.0, momentum_old, 0, dt*b1, momentum_k1, 0, 0, 2, 0);
446 200 : 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 200 : amrex::MultiFab::Saxpy(density_new, dt*b2, density_k2, 0, 0, 1, 0);
449 200 : amrex::MultiFab::Saxpy(momentum_new, dt*b2, momentum_k2, 0, 0, 2, 0);
450 200 : amrex::MultiFab::Saxpy(energy_new, dt*b2, energy_k2, 0, 0, 1, 0);
451 : // 2. ynew += dt*b3*k3
452 200 : amrex::MultiFab::Saxpy(density_new, dt*b3, density_k3, 0, 0, 1, 0);
453 200 : amrex::MultiFab::Saxpy(momentum_new, dt*b3, momentum_k3, 0, 0, 2, 0);
454 200 : amrex::MultiFab::Saxpy(energy_new, dt*b3, energy_k3, 0, 0, 1, 0);
455 :
456 200 : }
457 :
458 :
459 :
460 200 : else if (scheme == IntegrationScheme::RK4)
461 : {
462 : //
463 : // RK4 time integration scheme:
464 : //
465 :
466 200 : const amrex::BoxArray &ba = density_mf[lev]->boxArray();
467 200 : const amrex::DistributionMapping &dm = density_mf[lev]->DistributionMap();
468 200 : const int ng = density_mf[lev]->nGrow();
469 :
470 : // handles to old solution
471 200 : const amrex::MultiFab &density_old = *density_old_mf[lev];
472 200 : const amrex::MultiFab &momentum_old = *momentum_old_mf[lev];
473 200 : const amrex::MultiFab &energy_old = *energy_old_mf[lev];
474 :
475 : // runge kutta stages
476 200 : amrex::MultiFab density_k1(ba,dm,1,0), momentum_k1(ba,dm,2,0), energy_k1(ba,dm,1,0);
477 200 : amrex::MultiFab density_k2(ba,dm,1,0), momentum_k2(ba,dm,2,0), energy_k2(ba,dm,1,0);
478 200 : amrex::MultiFab density_k3(ba,dm,1,0), momentum_k3(ba,dm,2,0), energy_k3(ba,dm,1,0);
479 200 : 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 200 : 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 200 : density_st.ParallelCopyToGhost(*density_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
486 200 : momentum_st.ParallelCopyToGhost(*momentum_old_mf[lev],0,0,2,amrex::IntVect(1),amrex::IntVect(1));
487 200 : energy_st.ParallelCopyToGhost(*energy_old_mf[lev],0,0,1,amrex::IntVect(1),amrex::IntVect(1));
488 :
489 :
490 : // handles to new solution
491 200 : amrex::MultiFab &density_new = *density_mf[lev];
492 200 : amrex::MultiFab &momentum_new = *momentum_mf[lev];
493 200 : amrex::MultiFab &energy_new = *energy_mf[lev];
494 :
495 :
496 : //
497 : // K1
498 : //
499 :
500 200 : 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 200 : amrex::MultiFab::LinComb(density_st, 1.0, density_old, 0, dt/2.0, density_k1, 0, 0, 1, 0);
510 200 : amrex::MultiFab::LinComb(momentum_st, 1.0, momentum_old, 0, dt/2.0, momentum_k1, 0, 0, 2, 0);
511 200 : amrex::MultiFab::LinComb(energy_st, 1.0, energy_old, 0, dt/2.0, energy_k1, 0, 0, 1, 0);
512 :
513 200 : density_bc->FillBoundary(density_st, 0, 1, time, 0); density_st.FillBoundary(true);
514 200 : momentum_bc->FillBoundary(momentum_st, 0, 2, time, 0); momentum_st.FillBoundary(true);
515 200 : energy_bc->FillBoundary(energy_st,0,1,time,0); energy_st.FillBoundary(true);
516 :
517 :
518 200 : 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 200 : amrex::MultiFab::LinComb(density_st, 1.0, density_old, 0, dt/2.0, density_k2, 0, 0, 1, 0);
528 200 : amrex::MultiFab::LinComb(momentum_st, 1.0, momentum_old, 0, dt/2.0, momentum_k2, 0, 0, 2, 0);
529 200 : amrex::MultiFab::LinComb(energy_st, 1.0, energy_old, 0, dt/2.0, energy_k2, 0, 0, 1, 0);
530 :
531 200 : density_bc->FillBoundary(density_st, 0, 1, time, 0); density_st.FillBoundary(true);
532 200 : momentum_bc->FillBoundary(momentum_st, 0, 2, time, 0); momentum_st.FillBoundary(true);
533 200 : energy_bc->FillBoundary(energy_st,0,1,time,0); energy_st.FillBoundary(true);
534 :
535 200 : 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 200 : amrex::MultiFab::LinComb(density_st, 1.0, density_old, 0, dt, density_k3, 0, 0, 1, 0);
545 200 : amrex::MultiFab::LinComb(momentum_st, 1.0, momentum_old, 0, dt, momentum_k3, 0, 0, 2, 0);
546 200 : amrex::MultiFab::LinComb(energy_st, 1.0, energy_old, 0, dt, energy_k3, 0, 0, 1, 0);
547 :
548 200 : density_bc-> FillBoundary(density_st, 0, 1, time, 0); density_st.FillBoundary(true);
549 200 : momentum_bc->FillBoundary(momentum_st, 0, 2, time, 0); momentum_st.FillBoundary(true);
550 200 : energy_bc-> FillBoundary(energy_st, 0, 1, time, 0); energy_st.FillBoundary(true);
551 :
552 :
553 200 : 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 200 : amrex::MultiFab::LinComb(density_new, 1.0, density_old, 0, (dt/6.0), density_k1, 0, 0, 1, 0);
560 200 : amrex::MultiFab::LinComb(momentum_new, 1.0, momentum_old, 0, (dt/6.0), momentum_k1, 0, 0, 2, 0);
561 200 : 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 200 : amrex::MultiFab::Saxpy(density_new, (dt/3.0), density_k2, 0, 0, 1, 0);
565 200 : amrex::MultiFab::Saxpy(momentum_new, (dt/3.0), momentum_k2, 0, 0, 2, 0);
566 200 : amrex::MultiFab::Saxpy(energy_new, (dt/3.0), energy_k2, 0, 0, 1, 0);
567 :
568 : // [new] += (2 dt/6)k3
569 200 : amrex::MultiFab::Saxpy(density_new, (dt/3.0), density_k3, 0, 0, 1, 0);
570 200 : amrex::MultiFab::Saxpy(momentum_new, (dt/3.0), momentum_k3, 0, 0, 2, 0);
571 200 : amrex::MultiFab::Saxpy(energy_new, (dt/3.0), energy_k3, 0, 0, 1, 0);
572 :
573 : // [new] += (dt/6)k4
574 200 : amrex::MultiFab::Saxpy(density_new, (dt/6.0), density_k4, 0, 0, 1, 0);
575 200 : amrex::MultiFab::Saxpy(momentum_new, (dt/6.0), momentum_k4, 0, 0, 2, 0);
576 200 : amrex::MultiFab::Saxpy(energy_new, (dt/6.0), energy_k4, 0, 0, 1, 0);
577 :
578 200 : }
579 :
580 :
581 7300 : for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
582 : {
583 3650 : const amrex::Box& bx = mfi.validbox();
584 3650 : const Set::Scalar* DX = geom[lev].CellSize();
585 :
586 3650 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
587 3650 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
588 3650 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
589 3650 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
590 :
591 3650 : Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
592 3650 : Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
593 3650 : Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
594 :
595 3650 : Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
596 :
597 3650 : Set::Patch<Set::Scalar> u = velocity_mf.Patch(lev,mfi);
598 3650 : Set::Patch<Set::Scalar> Source = Source_mf.Patch(lev,mfi);
599 :
600 3650 : Set::Scalar *dt_max_handle = &dt_max;
601 :
602 3650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
603 : {
604 14182400 : if (eta(i,j,k) < cutoff)
605 : {
606 0 : rho_new(i,j,k,0) = rho_solid(i,j,k,0);
607 0 : M_new(i,j,k,0) = M_solid(i,j,k,0);
608 0 : M_new(i,j,k,1) = M_solid(i,j,k,1);
609 0 : E_new(i,j,k,0) = E_solid(i,j,k,0);
610 : }
611 :
612 7091200 : Set::Matrix gradu = Numeric::Gradient(u, i, j, k, DX);
613 14182400 : omega(i, j, k) = eta(i, j, k) * (gradu(1,0) - gradu(0,1));
614 :
615 7091200 : if (dynamictimestep.on)
616 : {
617 0 : *dt_max_handle = std::fabs(cfl * DX[0] / (u(i,j,k,0)*eta(i,j,k) + small));
618 0 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl * DX[1] / (u(i,j,k,1)*eta(i,j,k) + small)));
619 0 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[0]*DX[0] / (Source(i,j,k,1)+small)));
620 0 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[1]*DX[1] / (Source(i,j,k,2)+small)));
621 : }
622 7091200 : });
623 3650 : }
624 :
625 :
626 3650 : if (dynamictimestep.on)
627 : {
628 0 : this->DynamicTimestep_SyncTimeStep(lev,dt_max);
629 : }
630 :
631 3650 : }//end Advance
632 :
633 :
634 4650 : void 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 9300 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
644 : {
645 4650 : const amrex::Box& bx = mfi.growntilebox();
646 4650 : amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
647 :
648 4650 : Set::Patch<const Set::Scalar> rho = rho_mf.array(mfi); // density
649 4650 : Set::Patch<const Set::Scalar> M = M_mf.array(mfi); // momentum
650 4650 : 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 4650 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
653 4650 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
654 4650 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
655 :
656 4650 : Set::Patch<Set::Scalar> v = velocity_mf.Patch(lev,mfi);
657 4650 : Set::Patch<Set::Scalar> p = pressure_mf.Patch(lev,mfi);
658 :
659 4650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
660 : {
661 :
662 41203800 : Set::Scalar etarho_fluid = rho(i,j,k) - (1.-eta(i,j,k)) * rho_solid(i,j,k);
663 41203800 : Set::Scalar etaE_fluid = E(i,j,k) - (1.-eta(i,j,k)) * E_solid(i,j,k);
664 :
665 27469200 : Set::Vector etaM_fluid( M(i,j,k,0) - (1.-eta(i,j,k)) * M_solid(i,j,k,0),
666 68673000 : 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 13734600 : v(i,j,k,0) = etaM_fluid(0) / (etarho_fluid + small);
671 13734600 : v(i,j,k,1) = etaM_fluid(1) / (etarho_fluid + small);
672 41203800 : 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 13734600 : });
674 4650 : }
675 :
676 4650 : const Set::Scalar* DX = geom[lev].CellSize();
677 4650 : amrex::Box domain = geom[lev].Domain();
678 :
679 9300 : for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
680 : {
681 4650 : const amrex::Box& bx = mfi.validbox();
682 :
683 : // Inputs
684 4650 : Set::Patch<const Set::Scalar> rho = rho_mf.array(mfi);
685 4650 : Set::Patch<const Set::Scalar> E = E_mf.array(mfi);
686 4650 : Set::Patch<const Set::Scalar> M = M_mf.array(mfi);
687 :
688 : // Outputs
689 4650 : Set::Patch<Set::Scalar> rho_rhs = rho_rhs_mf.array(mfi);
690 4650 : Set::Patch<Set::Scalar> M_rhs = M_rhs_mf.array(mfi);
691 4650 : 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 4650 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
699 4650 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
700 4650 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
701 :
702 4650 : Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
703 :
704 4650 : Set::Patch<const Set::Scalar> eta = eta_old_mf.Patch(lev,mfi);
705 4650 : Set::Patch<const Set::Scalar> etadot = etadot_mf.Patch(lev,mfi);
706 4650 : Set::Patch<const Set::Scalar> velocity = velocity_mf.Patch(lev,mfi);
707 :
708 4650 : Set::Patch<const Set::Scalar> m0 = m0_mf.Patch(lev,mfi);
709 4650 : Set::Patch<const Set::Scalar> q = q_mf.Patch(lev,mfi);
710 4650 : Set::Patch<const Set::Scalar> _u0 = u0_mf.Patch(lev,mfi);
711 :
712 4650 : amrex::Array4<Set::Scalar> const& Source = (*Source_mf[lev]).array(mfi);
713 :
714 4650 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
715 : {
716 9139200 : auto sten = Numeric::GetStencil(i, j, k, domain);
717 :
718 : //Diffuse Sources
719 9139200 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
720 9139200 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
721 9139200 : Set::Matrix hess_eta = Numeric::Hessian(eta, i, j, k, 0, DX);
722 :
723 27417600 : Set::Vector u = Set::Vector(velocity(i, j, k, 0), velocity(i, j, k, 1));
724 27417600 : Set::Vector u0 = Set::Vector(_u0(i, j, k, 0), _u0(i, j, k, 1));
725 :
726 9139200 : Set::Matrix gradM = Numeric::Gradient(M, i, j, k, DX);
727 9139200 : Set::Vector gradrho = Numeric::Gradient(rho,i,j,k,0,DX);
728 9139200 : Set::Matrix hess_rho = Numeric::Hessian(rho,i,j,k,0,DX,sten);
729 18278400 : Set::Matrix gradu = (gradM - u*gradrho.transpose()) / rho(i,j,k);
730 :
731 27417600 : Set::Vector q0 = Set::Vector(q(i,j,k,0),q(i,j,k,1));
732 :
733 9139200 : Set::Scalar mdot0 = -m0(i,j,k) * grad_eta_mag;
734 9139200 : Set::Vector Pdot0 = Set::Vector::Zero();
735 9139200 : Set::Scalar qdot0 = q0.dot(grad_eta);
736 :
737 : // sten is necessary here because sometimes corner ghost
738 : // cells don't get filled
739 9139200 : Set::Matrix3 hess_M = Numeric::Hessian(M,i,j,k,DX);
740 9139200 : Set::Matrix3 hess_u = Set::Matrix3::Zero();
741 27417600 : for (int p = 0; p < 2; p++)
742 54835200 : for (int q = 0; q < 2; q++)
743 109670400 : for (int r = 0; r < 2; r++)
744 : {
745 73113600 : hess_u(r,p,q) =
746 73113600 : (hess_M(r,p,q) - gradu(r,q)*gradrho(p) - gradu(r,p)*gradrho(q) - u(r)*hess_rho(p,q))
747 146227200 : / rho(i,j,k);
748 : }
749 :
750 9139200 : Set::Vector Ldot0 = Set::Vector::Zero();
751 9139200 : Set::Vector div_tau = Set::Vector::Zero();
752 27417600 : for (int p = 0; p<2; p++)
753 54835200 : for (int q = 0; q<2; q++)
754 109670400 : for (int r = 0; r<2; r++)
755 219340800 : for (int s = 0; s<2; s++)
756 : {
757 146227200 : Set::Scalar Mpqrs = 0.0;
758 146227200 : if (p==r && q==s) Mpqrs += 0.5 * mu;
759 :
760 146227200 : Ldot0(p) += 0.5*Mpqrs * (u(r) - u0(r)) * hess_eta(q, s);
761 292454400 : div_tau(p) += 2.0*Mpqrs * hess_u(r,s,q);
762 : }
763 :
764 9139200 : Source(i,j, k, 0) = mdot0;
765 9139200 : Source(i,j, k, 1) = Pdot0(0) - Ldot0(0);
766 9139200 : Source(i,j, k, 2) = Pdot0(1) - Ldot0(1);
767 9139200 : 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 9139200 : Source(i,j,k,1) -= lagrange*(u-u0).dot(grad_eta)*grad_eta(0);
771 9139200 : Source(i,j,k,2) -= lagrange*(u-u0).dot(grad_eta)*grad_eta(1);
772 :
773 : //Godunov flux
774 : //states of total fields
775 9139200 : const int X = 0, Y = 1;
776 9139200 : Solver::Local::Riemann::State state_xlo(rho, M, E, i-1, j, k, X);
777 9139200 : Solver::Local::Riemann::State state_x (rho, M, E, i , j, k, X);
778 9139200 : Solver::Local::Riemann::State state_xhi(rho, M, E, i+1, j, k, X);
779 :
780 9139200 : Solver::Local::Riemann::State state_ylo(rho, M, E, i, j-1, k, Y);
781 9139200 : Solver::Local::Riemann::State state_y (rho, M, E, i, j , k, Y);
782 9139200 : Solver::Local::Riemann::State state_yhi(rho, M, E, i, j+1, k, Y);
783 :
784 : //states of solid fields
785 9139200 : Solver::Local::Riemann::State state_xlo_solid(rho_solid, M_solid, E_solid, i-1, j, k, X);
786 9139200 : Solver::Local::Riemann::State state_x_solid (rho_solid, M_solid, E_solid, i , j, k, X);
787 9139200 : Solver::Local::Riemann::State state_xhi_solid(rho_solid, M_solid, E_solid, i+1, j, k, X);
788 :
789 9139200 : Solver::Local::Riemann::State state_ylo_solid(rho_solid, M_solid, E_solid, i, j-1, k, Y);
790 9139200 : Solver::Local::Riemann::State state_y_solid (rho_solid, M_solid, E_solid, i, j , k, Y);
791 9139200 : Solver::Local::Riemann::State state_yhi_solid(rho_solid, M_solid, E_solid, i, j+1, k, Y);
792 :
793 27417600 : 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 27417600 : Solver::Local::Riemann::State state_x_fluid = (state_x - (1.0 - eta(i,j,k) )*state_x_solid ) / (eta(i,j,k) + small);
795 27417600 : 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 27417600 : 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 27417600 : Solver::Local::Riemann::State state_y_fluid = (state_y - (1.0 - eta(i,j,k) )*state_y_solid ) / (eta(i,j,k) + small);
799 27417600 : 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 9139200 : Solver::Local::Riemann::Flux flux_xlo, flux_ylo, flux_xhi, flux_yhi;
802 :
803 : try
804 : {
805 : //lo interface fluxes
806 18278400 : flux_xlo = riemannsolver->Solve(state_xlo_fluid, state_x_fluid, gamma, pref, small) * eta(i,j,k);
807 18278400 : flux_ylo = riemannsolver->Solve(state_ylo_fluid, state_y_fluid, gamma, pref, small) * eta(i,j,k);
808 :
809 : //hi interface fluxes
810 18278400 : flux_xhi = riemannsolver->Solve(state_x_fluid, state_xhi_fluid, gamma, pref, small) * eta(i,j,k);
811 18278400 : flux_yhi = riemannsolver->Solve(state_y_fluid, state_yhi_fluid, gamma, pref, small) * eta(i,j,k);
812 : }
813 0 : catch(...)
814 : {
815 0 : Util::ParallelMessage(INFO,"lev=",lev);
816 0 : Util::ParallelMessage(INFO,"i=",i,"j=",j);
817 0 : Util::Abort(INFO);
818 0 : }
819 :
820 :
821 : Set::Scalar drhof_dt =
822 9139200 : (flux_xlo.mass - flux_xhi.mass) / DX[0] +
823 9139200 : (flux_ylo.mass - flux_yhi.mass) / DX[1] +
824 9139200 : Source(i, j, k, 0);
825 :
826 18278400 : rho_rhs(i,j,k) =
827 : // rho_new(i, j, k) = rho(i, j, k) +
828 : //(
829 9139200 : drhof_dt +
830 : // todo add drhos_dt term if want time-evolving rhos
831 45696000 : 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 9139200 : (flux_xlo.momentum_normal - flux_xhi.momentum_normal ) / DX[0] +
839 18278400 : (flux_ylo.momentum_tangent - flux_yhi.momentum_tangent) / DX[1] +
840 9139200 : div_tau(0) * eta(i,j,k) +
841 9139200 : Source(i, j, k, 1);
842 :
843 18278400 : M_rhs(i,j,k,0) =
844 : //M_new(i, j, k, 0) = M(i, j, k, 0) +
845 : // (
846 9139200 : dMxf_dt +
847 : // todo add dMs_dt term if want time-evolving Ms
848 45696000 : 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 9139200 : (flux_xlo.momentum_tangent - flux_xhi.momentum_tangent) / DX[0] +
854 18278400 : (flux_ylo.momentum_normal - flux_yhi.momentum_normal ) / DX[1] +
855 9139200 : div_tau(1) * eta(i,j,k) +
856 9139200 : Source(i, j, k, 2);
857 :
858 18278400 : M_rhs(i,j,k,1) =
859 : //M_new(i, j, k, 1) = M(i, j, k, 1) +
860 : //(
861 9139200 : dMyf_dt +
862 : // todo add dMs_dt term if want time-evolving Ms
863 45696000 : 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 9139200 : (flux_xlo.energy - flux_xhi.energy) / DX[0] +
869 9139200 : (flux_ylo.energy - flux_yhi.energy) / DX[1] +
870 9139200 : Source(i, j, k, 3);
871 :
872 18278400 : E_rhs(i,j,k) =
873 : // E_new(i, j, k) = E(i, j, k) +
874 : // (
875 9139200 : dEf_dt +
876 : // todo add dEs_dt term if want time-evolving Es
877 45696000 : 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 18278400 : omega(i, j, k) = eta(i, j, k) * (gradu(1,0) - gradu(0,1));
945 9139200 : });
946 4650 : }
947 4650 : }
948 :
949 0 : void Hydro::Regrid(int lev, Set::Scalar /* time */)
950 : {
951 : BL_PROFILE("Integrator::Hydro::Regrid");
952 0 : Source_mf[lev]->setVal(0.0);
953 0 : if (lev < finest_level) return;
954 :
955 0 : 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)
959 0 : void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar, int)
960 : {
961 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
962 :
963 0 : const Set::Scalar* DX = geom[lev].CellSize();
964 0 : 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 0 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi) {
968 0 : const amrex::Box& bx = mfi.tilebox();
969 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
970 0 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
971 :
972 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
973 0 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
974 0 : if (grad_eta.lpNorm<2>() * dr * 2 > eta_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
975 0 : });
976 0 : }
977 :
978 : // Vorticity criterion for refinement
979 0 : for (amrex::MFIter mfi(*vorticity_mf[lev], true); mfi.isValid(); ++mfi) {
980 0 : const amrex::Box& bx = mfi.tilebox();
981 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
982 0 : amrex::Array4<const Set::Scalar> const& omega = (*vorticity_mf[lev]).array(mfi);
983 :
984 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
985 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
986 0 : Set::Vector grad_omega = Numeric::Gradient(omega, i, j, k, 0, DX, sten);
987 0 : if (grad_omega.lpNorm<2>() * dr * 2 > omega_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
988 0 : });
989 0 : }
990 :
991 : // Gradu criterion for refinement
992 0 : for (amrex::MFIter mfi(*velocity_mf[lev], true); mfi.isValid(); ++mfi) {
993 0 : const amrex::Box& bx = mfi.tilebox();
994 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
995 0 : amrex::Array4<const Set::Scalar> const& v = (*velocity_mf[lev]).array(mfi);
996 :
997 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
998 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
999 0 : Set::Matrix grad_u = Numeric::Gradient(v, i, j, k, DX, sten);
1000 0 : if (grad_u.lpNorm<2>() * dr * 2 > gradu_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
1001 0 : });
1002 0 : }
1003 :
1004 : // Pressure criterion for refinement
1005 0 : for (amrex::MFIter mfi(*pressure_mf[lev], true); mfi.isValid(); ++mfi) {
1006 0 : const amrex::Box& bx = mfi.tilebox();
1007 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
1008 0 : amrex::Array4<const Set::Scalar> const& p = (*pressure_mf[lev]).array(mfi);
1009 :
1010 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
1011 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
1012 0 : Set::Vector grad_p = Numeric::Gradient(p, i, j, k, 0, DX, sten);
1013 0 : if (grad_p.lpNorm<2>() * dr * 2 > p_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
1014 0 : });
1015 0 : }
1016 :
1017 : // Density criterion for refinement
1018 0 : for (amrex::MFIter mfi(*density_mf[lev], true); mfi.isValid(); ++mfi) {
1019 0 : const amrex::Box& bx = mfi.tilebox();
1020 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
1021 0 : amrex::Array4<const Set::Scalar> const& rho = (*density_mf[lev]).array(mfi);
1022 :
1023 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
1024 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
1025 0 : Set::Vector grad_rho = Numeric::Gradient(rho, i, j, k, 0, DX, sten);
1026 0 : if (grad_rho.lpNorm<2>() * dr * 2 > rho_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
1027 0 : });
1028 0 : }
1029 :
1030 0 : }
1031 :
1032 : }
1033 :
1034 :
1035 : #endif
|
