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 :
14 : #if AMREX_SPACEDIM == 2
15 :
16 : namespace Integrator
17 : {
18 :
19 4 : Hydro::Hydro(IO::ParmParse& pp) : Hydro()
20 : {
21 12 : pp.queryclass(*this);
22 4 : }
23 :
24 : void
25 4 : Hydro::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 : // eta-based refinement
36 24 : pp.query_default("eta_refinement_criterion", value.eta_refinement_criterion , 0.01);
37 : // vorticity-based refinement
38 24 : pp.query_default("omega_refinement_criterion", value.omega_refinement_criterion, 0.01);
39 : // velocity gradient-based refinement
40 24 : pp.query_default("gradu_refinement_criterion", value.gradu_refinement_criterion, 0.01);
41 : // pressure-based refinement
42 24 : pp.query_default("p_refinement_criterion", value.p_refinement_criterion, 1e100);
43 : // density-based refinement
44 24 : pp.query_default("rho_refinement_criterion", value.rho_refinement_criterion, 1e100);
45 :
46 24 : pp_query_required("gamma", value.gamma); // gamma for gamma law
47 24 : pp_query_required("cfl", value.cfl); // cfl condition
48 24 : pp_query_default("cfl_v", value.cfl_v,1E100); // cfl condition
49 24 : pp_query_required("mu", value.mu); // linear viscosity coefficient
50 32 : pp_forbid("Lfactor","replaced with mu");
51 : //pp_query_default("Lfactor", value.Lfactor,1.0); // (to be removed) test factor for viscous source
52 32 : pp_forbid("Pfactor","replaced with mu");
53 : //pp_query_default("Pfactor", value.Pfactor,1.0); // (to be removed) test factor for viscous source
54 24 : pp_query_default("pref", value.pref,1.0); // reference pressure for Roe solver
55 :
56 32 : pp_forbid("rho.bc","--> density.bc");
57 32 : pp_forbid("p.bc","--> pressure.bc");
58 32 : pp_forbid("v.bc", "--> velocity.bc");
59 32 : pp_forbid("pressure.bc","--> energy.bc");
60 28 : pp_forbid("velocity.bc","--> momentum.bc");
61 :
62 : // Boundary condition for density
63 8 : pp.select_default<BC::Constant,BC::Expression>("density.bc",value.density_bc,1);
64 : // Boundary condition for energy
65 8 : pp.select_default<BC::Constant,BC::Expression>("energy.bc",value.energy_bc,1);
66 : // Boundary condition for momentum
67 8 : pp.select_default<BC::Constant,BC::Expression>("momentum.bc",value.momentum_bc,2);
68 : // Boundary condition for phase field order parameter
69 12 : pp.select_default<BC::Constant,BC::Expression>("pf.eta.bc",value.eta_bc,1);
70 :
71 24 : pp_query_default("small",value.small,1E-8); // small regularization value
72 24 : pp_query_default("cutoff",value.cutoff,-1E100); // cutoff value
73 24 : pp_query_default("lagrange",value.lagrange,0.0); // lagrange no-penetration factor
74 :
75 28 : pp_forbid("roefix","--> solver.roe.entropy_fix"); // Roe solver entropy fix
76 :
77 : }
78 : // Register FabFields:
79 : {
80 4 : int nghost = 2;
81 :
82 12 : value.RegisterNewFab(value.eta_mf, value.eta_bc, 1, nghost, "eta", true );
83 12 : value.RegisterNewFab(value.eta_old_mf, value.eta_bc, 1, nghost, "eta_old", true);
84 12 : value.RegisterNewFab(value.etadot_mf, value.eta_bc, 1, nghost, "etadot", true );
85 :
86 12 : value.RegisterNewFab(value.density_mf, value.density_bc, 1, nghost, "density", true );
87 12 : value.RegisterNewFab(value.density_old_mf, value.density_bc, 1, nghost, "density_old", false);
88 :
89 12 : value.RegisterNewFab(value.energy_mf, value.energy_bc, 1, nghost, "energy", true );
90 12 : value.RegisterNewFab(value.energy_old_mf, value.energy_bc, 1, nghost, "energy_old" , false);
91 :
92 16 : value.RegisterNewFab(value.momentum_mf, value.momentum_bc, 2, nghost, "momentum", true , {"x","y"});
93 12 : value.RegisterNewFab(value.momentum_old_mf, value.momentum_bc, 2, nghost, "momentum_old", false);
94 :
95 12 : value.RegisterNewFab(value.pressure_mf, &value.bc_nothing, 1, nghost, "pressure", true);
96 16 : value.RegisterNewFab(value.velocity_mf, &value.bc_nothing, 2, nghost, "velocity", true, {"x","y"});
97 12 : value.RegisterNewFab(value.vorticity_mf, &value.bc_nothing, 1, nghost, "vorticity", true);
98 :
99 12 : value.RegisterNewFab(value.m0_mf, &value.bc_nothing, 1, 0, "m0", true);
100 16 : value.RegisterNewFab(value.u0_mf, &value.bc_nothing, 2, 0, "u0", true, {"x","y"});
101 16 : value.RegisterNewFab(value.q_mf, &value.bc_nothing, 2, 0, "q", true, {"x","y"});
102 :
103 16 : value.RegisterNewFab(value.solid.momentum_mf, &value.neumann_bc_D, 2, nghost, "solid.momentum", true, {"x","y"});
104 12 : value.RegisterNewFab(value.solid.density_mf, &value.neumann_bc_1, 1, nghost, "solid.density", true);
105 12 : value.RegisterNewFab(value.solid.energy_mf, &value.neumann_bc_1, 1, nghost, "solid.energy", true);
106 :
107 12 : value.RegisterNewFab(value.Source_mf, &value.bc_nothing, 4, 0, "Source", true);
108 : }
109 :
110 32 : pp_forbid("Velocity.ic.type", "--> velocity.ic.type");
111 32 : pp_forbid("Pressure.ic", "--> pressure.ic");
112 32 : pp_forbid("SolidMomentum.ic", "--> solid.momentum.ic");
113 32 : pp_forbid("SolidDensity.ic.type", "--> solid.density.ic.type");
114 32 : pp_forbid("SolidEnergy.ic.type", "--> solid.energy.ic.type");
115 32 : pp_forbid("Density.ic.type", "--> density.ic.type");
116 32 : pp_forbid("rho_injected.ic.type","no longer using rho_injected use m0 instead");
117 28 : pp.forbid("mdot.ic.type", "replace mdot with u0");
118 :
119 :
120 : // ORDER PARAMETER
121 :
122 : // eta initial condition
123 8 : pp.select_default<IC::Constant,IC::Laminate,IC::Expression,IC::BMP,IC::PNG>("eta.ic",value.eta_ic,value.geom);
124 :
125 : // PRIMITIVE FIELD INITIAL CONDITIONS
126 :
127 : // velocity initial condition
128 8 : pp.select_default<IC::Constant,IC::Expression>("velocity.ic",value.velocity_ic,value.geom);
129 : // solid pressure initial condition
130 8 : pp.select_default<IC::Constant,IC::Expression>("pressure.ic",value.pressure_ic,value.geom);
131 : // density initial condition type
132 8 : pp.select_default<IC::Constant,IC::Expression>("density.ic",value.density_ic,value.geom);
133 :
134 :
135 : // SOLID FIELDS
136 :
137 : // solid momentum initial condition
138 8 : pp.select_default<IC::Constant,IC::Expression>("solid.momentum.ic",value.solid.momentum_ic,value.geom);
139 : // solid density initial condition
140 8 : pp.select_default<IC::Constant,IC::Expression>("solid.density.ic",value.solid.density_ic,value.geom);
141 : // solid energy initial condition
142 8 : pp.select_default<IC::Constant,IC::Expression>("solid.energy.ic",value.solid.energy_ic,value.geom);
143 :
144 :
145 : // DIFFUSE BOUNDARY SOURCES
146 :
147 : // diffuse boundary prescribed mass flux
148 8 : pp.select_default<IC::Constant,IC::Expression>("m0.ic",value.ic_m0,value.geom);
149 : // diffuse boundary prescribed velocity
150 8 : pp.select_default<IC::Constant,IC::Expression>("u0.ic",value.ic_u0,value.geom);
151 : // diffuse boundary prescribed heat flux
152 8 : pp.select_default<IC::Constant,IC::Expression>("q.ic",value.ic_q,value.geom);
153 :
154 : // Riemann solver
155 8 : pp.select_default<Solver::Local::Riemann::Roe>("solver",value.roesolver);
156 4 : }
157 :
158 :
159 4 : void Hydro::Initialize(int lev)
160 : {
161 : BL_PROFILE("Integrator::Hydro::Initialize");
162 :
163 4 : eta_ic ->Initialize(lev, eta_mf, 0.0);
164 4 : eta_ic ->Initialize(lev, eta_old_mf, 0.0);
165 4 : etadot_mf[lev] ->setVal(0.0);
166 :
167 : //flux_mf[lev] ->setVal(0.0);
168 :
169 4 : velocity_ic ->Initialize(lev, velocity_mf, 0.0);
170 4 : pressure_ic ->Initialize(lev, pressure_mf, 0.0);
171 4 : density_ic ->Initialize(lev, density_mf, 0.0);
172 :
173 4 : density_ic ->Initialize(lev, density_old_mf, 0.0);
174 :
175 :
176 4 : solid.density_ic ->Initialize(lev, solid.density_mf, 0.0);
177 4 : solid.momentum_ic->Initialize(lev, solid.momentum_mf, 0.0);
178 4 : solid.energy_ic ->Initialize(lev, solid.energy_mf, 0.0);
179 :
180 4 : ic_m0 ->Initialize(lev, m0_mf, 0.0);
181 4 : ic_u0 ->Initialize(lev, u0_mf, 0.0);
182 4 : ic_q ->Initialize(lev, q_mf, 0.0);
183 :
184 4 : Source_mf[lev] ->setVal(0.0);
185 :
186 4 : Mix(lev);
187 4 : }
188 :
189 4 : void Hydro::Mix(int lev)
190 : {
191 8 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
192 : {
193 4 : const amrex::Box& bx = mfi.growntilebox();
194 :
195 4 : Set::Patch<const Set::Scalar> eta = eta_mf.Patch(lev,mfi);
196 :
197 4 : Set::Patch<const Set::Scalar> v = velocity_mf.Patch(lev,mfi);
198 4 : Set::Patch<const Set::Scalar> p = pressure_mf.Patch(lev,mfi);
199 4 : Set::Patch<Set::Scalar> rho = density_mf.Patch(lev,mfi);
200 4 : Set::Patch<Set::Scalar> rho_old = density_old_mf.Patch(lev,mfi);
201 4 : Set::Patch<Set::Scalar> M = momentum_mf.Patch(lev,mfi);
202 4 : Set::Patch<Set::Scalar> M_old = momentum_old_mf.Patch(lev,mfi);
203 4 : Set::Patch<Set::Scalar> E = energy_mf.Patch(lev,mfi);
204 4 : Set::Patch<Set::Scalar> E_old = energy_old_mf.Patch(lev,mfi);
205 4 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
206 4 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
207 4 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
208 :
209 :
210 4 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
211 : {
212 66000 : rho(i, j, k) = eta(i, j, k) * rho(i, j, k) + (1.0 - eta(i, j, k)) * rho_solid(i, j, k);
213 26400 : rho_old(i, j, k) = rho(i, j, k);
214 :
215 79200 : 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));
216 79200 : 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));
217 26400 : M_old(i, j, k, 0) = M(i, j, k, 0);
218 26400 : M_old(i, j, k, 1) = M(i, j, k, 1);
219 :
220 26400 : E(i, j, k) =
221 92400 : (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)
222 13200 : +
223 39600 : E_solid(i, j, k) * (1.0 - eta(i, j, k));
224 26400 : E_old(i, j, k) = E(i, j, k);
225 13200 : });
226 4 : }
227 4 : c_max = 0.0;
228 4 : vx_max = 0.0;
229 4 : vy_max = 0.0;
230 4 : }
231 :
232 3250 : void Hydro::UpdateEta(int lev, Set::Scalar time)
233 : {
234 3250 : eta_ic->Initialize(lev, eta_mf, time);
235 3250 : }
236 :
237 3250 : void Hydro::TimeStepBegin(Set::Scalar, int /*iter*/)
238 : {
239 :
240 3250 : }
241 :
242 3250 : void Hydro::TimeStepComplete(Set::Scalar, int lev)
243 : {
244 3250 : Integrator::DynamicTimestep_Update();
245 :
246 3250 : return;
247 :
248 : const Set::Scalar* DX = geom[lev].CellSize();
249 :
250 : amrex::ParallelDescriptor::ReduceRealMax(c_max);
251 : amrex::ParallelDescriptor::ReduceRealMax(vx_max);
252 : amrex::ParallelDescriptor::ReduceRealMax(vy_max);
253 :
254 : Set::Scalar new_timestep = cfl / ((c_max + vx_max) / DX[0] + (c_max + vy_max) / DX[1]);
255 :
256 : Util::Assert(INFO, TEST(AMREX_SPACEDIM == 2));
257 :
258 : SetTimestep(new_timestep);
259 : }
260 :
261 3250 : void Hydro::Advance(int lev, Set::Scalar time, Set::Scalar dt)
262 : {
263 :
264 3250 : std::swap(eta_old_mf, eta_mf);
265 3250 : std::swap(density_old_mf[lev], density_mf[lev]);
266 3250 : std::swap(momentum_old_mf[lev], momentum_mf[lev]);
267 3250 : std::swap(energy_old_mf[lev], energy_mf[lev]);
268 3250 : Set::Scalar dt_max = std::numeric_limits<Set::Scalar>::max();
269 :
270 3250 : UpdateEta(lev, time);
271 :
272 6500 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
273 : {
274 3250 : const amrex::Box& bx = mfi.growntilebox();
275 3250 : amrex::Array4<const Set::Scalar> const& eta_new = (*eta_mf[lev]).array(mfi);
276 3250 : amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
277 3250 : amrex::Array4<Set::Scalar> const& etadot = (*etadot_mf[lev]).array(mfi);
278 :
279 3250 : Set::Patch<const Set::Scalar> rho = density_old_mf.Patch(lev,mfi);
280 3250 : Set::Patch<const Set::Scalar> M = momentum_old_mf.Patch(lev,mfi);
281 3250 : Set::Patch<const Set::Scalar> E = energy_old_mf.Patch(lev,mfi); // total energy (internal energy + kinetic energy) per unit volume
282 : // E/rho = e + 0.5*v^2
283 :
284 3250 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
285 3250 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
286 3250 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
287 :
288 3250 : Set::Patch<Set::Scalar> v = velocity_mf.Patch(lev,mfi);
289 3250 : Set::Patch<Set::Scalar> p = pressure_mf.Patch(lev,mfi);
290 :
291 3250 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
292 : {
293 37764000 : etadot(i, j, k) = (eta_new(i, j, k) - eta(i, j, k)) / dt;
294 :
295 37764000 : Set::Scalar etarho_fluid = rho(i,j,k) - (1.-eta(i,j,k)) * rho_solid(i,j,k);
296 37764000 : Set::Scalar etaE_fluid = E(i,j,k) - (1.-eta(i,j,k)) * E_solid(i,j,k);
297 :
298 25176000 : Set::Vector etaM_fluid( M(i,j,k,0) - (1.-eta(i,j,k)) * M_solid(i,j,k,0),
299 62940000 : M(i,j,k,1) - (1.-eta(i,j,k)) * M_solid(i,j,k,1) );
300 :
301 : //THESE ARE FLUID VELOCITY AND PRESSURE
302 :
303 12588000 : v(i,j,k,0) = etaM_fluid(0) / (etarho_fluid + small);
304 12588000 : v(i,j,k,1) = etaM_fluid(1) / (etarho_fluid + small);
305 :
306 37764000 : 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;
307 12588000 : });
308 3250 : }
309 :
310 3250 : const Set::Scalar* DX = geom[lev].CellSize();
311 3250 : amrex::Box domain = geom[lev].Domain();
312 :
313 6500 : for (amrex::MFIter mfi(*eta_mf[lev], false); mfi.isValid(); ++mfi)
314 : {
315 3250 : const amrex::Box& bx = mfi.validbox();
316 :
317 3250 : Set::Patch<const Set::Scalar> rho = density_old_mf.Patch(lev,mfi);
318 3250 : Set::Patch<const Set::Scalar> E = energy_old_mf.Patch(lev,mfi);
319 3250 : Set::Patch<const Set::Scalar> M = momentum_old_mf.Patch(lev,mfi);
320 :
321 3250 : Set::Patch<Set::Scalar> rho_new = density_mf.Patch(lev,mfi);
322 3250 : Set::Patch<Set::Scalar> E_new = energy_mf.Patch(lev,mfi);
323 3250 : Set::Patch<Set::Scalar> M_new = momentum_mf.Patch(lev,mfi);
324 :
325 3250 : Set::Patch<const Set::Scalar> rho_solid = solid.density_mf.Patch(lev,mfi);
326 3250 : Set::Patch<const Set::Scalar> M_solid = solid.momentum_mf.Patch(lev,mfi);
327 3250 : Set::Patch<const Set::Scalar> E_solid = solid.energy_mf.Patch(lev,mfi);
328 :
329 3250 : Set::Patch<Set::Scalar> omega = vorticity_mf.Patch(lev,mfi);
330 :
331 3250 : Set::Patch<const Set::Scalar> eta = eta_old_mf.Patch(lev,mfi);
332 3250 : Set::Patch<const Set::Scalar> etadot = etadot_mf.Patch(lev,mfi);
333 3250 : Set::Patch<const Set::Scalar> velocity = velocity_mf.Patch(lev,mfi);
334 :
335 3250 : Set::Patch<const Set::Scalar> m0 = m0_mf.Patch(lev,mfi);
336 3250 : Set::Patch<const Set::Scalar> q = q_mf.Patch(lev,mfi);
337 3250 : Set::Patch<const Set::Scalar> _u0 = u0_mf.Patch(lev,mfi);
338 :
339 3250 : amrex::Array4<Set::Scalar> const& Source = (*Source_mf[lev]).array(mfi);
340 :
341 3250 : Set::Scalar *dt_max_handle = &dt_max;
342 :
343 3250 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
344 : {
345 6272000 : auto sten = Numeric::GetStencil(i, j, k, domain);
346 :
347 : //Diffuse Sources
348 6272000 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
349 6272000 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
350 6272000 : Set::Matrix hess_eta = Numeric::Hessian(eta, i, j, k, 0, DX);
351 :
352 18816000 : Set::Vector u = Set::Vector(velocity(i, j, k, 0), velocity(i, j, k, 1));
353 18816000 : Set::Vector u0 = Set::Vector(_u0(i, j, k, 0), _u0(i, j, k, 1));
354 :
355 6272000 : Set::Matrix gradM = Numeric::Gradient(M, i, j, k, DX);
356 6272000 : Set::Vector gradrho = Numeric::Gradient(rho,i,j,k,0,DX);
357 6272000 : Set::Matrix hess_rho = Numeric::Hessian(rho,i,j,k,0,DX,sten);
358 12544000 : Set::Matrix gradu = (gradM - u*gradrho.transpose()) / rho(i,j,k);
359 :
360 18816000 : Set::Vector q0 = Set::Vector(q(i,j,k,0),q(i,j,k,1));
361 :
362 :
363 6272000 : Set::Scalar mdot0 = -m0(i,j,k) * grad_eta_mag;
364 6272000 : Set::Vector Pdot0 = Set::Vector::Zero();
365 6272000 : Set::Scalar qdot0 = q0.dot(grad_eta);
366 :
367 6272000 : Set::Matrix3 hess_M = Numeric::Hessian(M,i,j,k,DX);
368 6272000 : Set::Matrix3 hess_u = Set::Matrix3::Zero();
369 18816000 : for (int p = 0; p < 2; p++)
370 37632000 : for (int q = 0; q < 2; q++)
371 75264000 : for (int r = 0; r < 2; r++)
372 : {
373 50176000 : hess_u(r,p,q) =
374 50176000 : (hess_M(r,p,q) - gradu(r,q)*gradrho(p) - gradu(r,p)*gradrho(q) - u(r)*hess_rho(p,q))
375 100352000 : / rho(i,j,k);
376 : }
377 :
378 6272000 : Set::Vector Ldot0 = Set::Vector::Zero();
379 6272000 : Set::Vector div_tau = Set::Vector::Zero();
380 18816000 : for (int p = 0; p<2; p++)
381 37632000 : for (int q = 0; q<2; q++)
382 75264000 : for (int r = 0; r<2; r++)
383 150528000 : for (int s = 0; s<2; s++)
384 : {
385 100352000 : Set::Scalar Mpqrs = 0.0;
386 100352000 : if (p==r && q==s) Mpqrs += 0.5 * mu;
387 :
388 100352000 : Ldot0(p) += 0.5*Mpqrs * (u(r) - u0(r)) * hess_eta(q, s);
389 200704000 : div_tau(p) += 2.0*Mpqrs * hess_u(r,s,q);
390 : }
391 :
392 6272000 : Source(i,j, k, 0) = mdot0;
393 6272000 : Source(i,j, k, 1) = Pdot0(0) - Ldot0(0);
394 6272000 : Source(i,j, k, 2) = Pdot0(1) - Ldot0(1);
395 6272000 : Source(i,j, k, 3) = qdot0;// - Ldot0(0)*v(i,j,k,0) - Ldot0(1)*v(i,j,k,1);
396 :
397 : // Lagrange terms to enforce no-penetration
398 6272000 : Source(i,j,k,1) -= lagrange*u.dot(grad_eta)*grad_eta(0);
399 6272000 : Source(i,j,k,2) -= lagrange*u.dot(grad_eta)*grad_eta(1);
400 :
401 : //Godunov flux
402 : //states of total fields
403 6272000 : const int X = 0, Y = 1;
404 6272000 : Solver::Local::Riemann::State state_xlo(rho, M, E, i-1, j, k, X);
405 6272000 : Solver::Local::Riemann::State state_x (rho, M, E, i , j, k, X);
406 6272000 : Solver::Local::Riemann::State state_xhi(rho, M, E, i+1, j, k, X);
407 :
408 6272000 : Solver::Local::Riemann::State state_ylo(rho, M, E, i, j-1, k, Y);
409 6272000 : Solver::Local::Riemann::State state_y (rho, M, E, i, j , k, Y);
410 6272000 : Solver::Local::Riemann::State state_yhi(rho, M, E, i, j+1, k, Y);
411 :
412 : //states of solid fields
413 6272000 : Solver::Local::Riemann::State state_xlo_solid(rho_solid, M_solid, E_solid, i-1, j, k, X);
414 6272000 : Solver::Local::Riemann::State state_x_solid (rho_solid, M_solid, E_solid, i , j, k, X);
415 6272000 : Solver::Local::Riemann::State state_xhi_solid(rho_solid, M_solid, E_solid, i+1, j, k, X);
416 :
417 6272000 : Solver::Local::Riemann::State state_ylo_solid(rho_solid, M_solid, E_solid, i, j-1, k, Y);
418 6272000 : Solver::Local::Riemann::State state_y_solid (rho_solid, M_solid, E_solid, i, j , k, Y);
419 6272000 : Solver::Local::Riemann::State state_yhi_solid(rho_solid, M_solid, E_solid, i, j+1, k, Y);
420 :
421 18816000 : 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);
422 18816000 : Solver::Local::Riemann::State state_x_fluid = (state_x - (1.0 - eta(i,j,k) )*state_x_solid ) / (eta(i,j,k) + small);
423 18816000 : 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);
424 :
425 18816000 : 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);
426 18816000 : Solver::Local::Riemann::State state_y_fluid = (state_y - (1.0 - eta(i,j,k) )*state_y_solid ) / (eta(i,j,k) + small);
427 18816000 : 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);
428 :
429 6272000 : Solver::Local::Riemann::Flux flux_xlo, flux_ylo, flux_xhi, flux_yhi;
430 :
431 : try
432 : {
433 : //lo interface fluxes
434 12544000 : flux_xlo = roesolver->Solve(state_xlo_fluid, state_x_fluid, gamma, pref, small) * eta(i,j,k);
435 12544000 : flux_ylo = roesolver->Solve(state_ylo_fluid, state_y_fluid, gamma, pref, small) * eta(i,j,k);
436 :
437 : //hi interface fluxes
438 12544000 : flux_xhi = roesolver->Solve(state_x_fluid, state_xhi_fluid, gamma, pref, small) * eta(i,j,k);
439 12544000 : flux_yhi = roesolver->Solve(state_y_fluid, state_yhi_fluid, gamma, pref, small) * eta(i,j,k);
440 : }
441 0 : catch(...)
442 : {
443 0 : Util::ParallelMessage(INFO,"lev=",lev);
444 0 : Util::ParallelMessage(INFO,"i=",i,"j=",j);
445 0 : Util::Abort(INFO);
446 0 : }
447 :
448 :
449 : Set::Scalar drhof_dt =
450 6272000 : (flux_xlo.mass - flux_xhi.mass) / DX[0] +
451 12544000 : (flux_ylo.mass - flux_yhi.mass) / DX[1] +
452 6272000 : Source(i, j, k, 0);
453 :
454 6272000 : rho_new(i, j, k) = rho(i, j, k) +
455 : (
456 6272000 : drhof_dt +
457 : // todo add drhos_dt term if want time-evolving rhos
458 25088000 : etadot(i,j,k) * (rho(i,j,k) - rho_solid(i,j,k)) / (eta(i,j,k) + small)
459 6272000 : ) * dt;
460 :
461 18816000 : if (rho_new(i,j,k) != rho_new(i,j,k))
462 : {
463 0 : Util::ParallelMessage(INFO,"lev=",lev);
464 0 : Util::ParallelMessage(INFO,"i=",i,"j=",j);
465 0 : Util::ParallelMessage(INFO,"drhof_dt",drhof_dt); // dies
466 0 : Util::ParallelMessage(INFO,"flux_xlo.mass",flux_xlo.mass);
467 0 : 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
468 0 : Util::ParallelMessage(INFO,"flux_ylo.mass",flux_ylo.mass);
469 0 : Util::ParallelMessage(INFO,"flux_xhi.mass",flux_yhi.mass);
470 0 : Util::ParallelMessage(INFO,"eta",eta(i,j,k));
471 0 : Util::ParallelMessage(INFO,"Source",Source(i,j,k,0));
472 0 : Util::ParallelMessage(INFO,"state_x",state_x); // <<<<
473 0 : Util::ParallelMessage(INFO,"state_y",state_y);
474 0 : Util::ParallelMessage(INFO,"state_x_solid",state_x_solid); // <<<<
475 0 : Util::ParallelMessage(INFO,"state_y_solid",state_y_solid);
476 0 : Util::ParallelMessage(INFO,"state_xhi",state_xhi); // <<<<
477 0 : Util::ParallelMessage(INFO,"state_yhi",state_yhi);
478 0 : Util::ParallelMessage(INFO,"state_xhi_solid",state_xhi_solid);
479 0 : Util::ParallelMessage(INFO,"state_yhi_solids",state_yhi_solid);
480 0 : Util::ParallelMessage(INFO,"state_xlo",state_xlo);
481 0 : Util::ParallelMessage(INFO,"state_ylo",state_ylo);
482 0 : Util::ParallelMessage(INFO,"state_xlo_solid",state_xlo_solid);
483 0 : Util::ParallelMessage(INFO,"state_ylo_solid",state_ylo_solid);
484 0 : Util::Exception(INFO);
485 : }
486 :
487 :
488 : Set::Scalar dMxf_dt =
489 6272000 : (flux_xlo.momentum_normal - flux_xhi.momentum_normal ) / DX[0] +
490 12544000 : (flux_ylo.momentum_tangent - flux_yhi.momentum_tangent) / DX[1] +
491 6272000 : div_tau(0) * eta(i,j,k) +
492 : //(mu * (lap_ux * eta(i, j, k))) +
493 6272000 : Source(i, j, k, 1);
494 :
495 6272000 : M_new(i, j, k, 0) = M(i, j, k, 0) +
496 : (
497 6272000 : dMxf_dt +
498 : // todo add dMs_dt term if want time-evolving Ms
499 25088000 : etadot(i,j,k)*(M(i,j,k,0) - M_solid(i,j,k,0)) / (eta(i,j,k) + small)
500 6272000 : ) * dt;
501 :
502 : Set::Scalar dMyf_dt =
503 6272000 : (flux_xlo.momentum_tangent - flux_xhi.momentum_tangent) / DX[0] +
504 12544000 : (flux_ylo.momentum_normal - flux_yhi.momentum_normal ) / DX[1] +
505 6272000 : div_tau(1) * eta(i,j,k) +
506 : //(mu * (lap_uy * eta(i, j, k))) +
507 6272000 : Source(i, j, k, 2);
508 :
509 6272000 : M_new(i, j, k, 1) = M(i, j, k, 1) +
510 : (
511 6272000 : dMyf_dt +
512 : // todo add dMs_dt term if want time-evolving Ms
513 25088000 : etadot(i,j,k)*(M(i,j,k,1) - M_solid(i,j,k,1)) / (eta(i,j,k)+small)
514 6272000 : )*dt;
515 :
516 : Set::Scalar dEf_dt =
517 6272000 : (flux_xlo.energy - flux_xhi.energy) / DX[0] +
518 6272000 : (flux_ylo.energy - flux_yhi.energy) / DX[1] +
519 6272000 : Source(i, j, k, 3);
520 :
521 6272000 : E_new(i, j, k) = E(i, j, k) +
522 : (
523 6272000 : dEf_dt +
524 : // todo add dEs_dt term if want time-evolving Es
525 25088000 : etadot(i,j,k)*(E(i,j,k) - E_solid(i,j,k)) / (eta(i,j,k)+small)
526 6272000 : ) * dt;
527 :
528 :
529 12544000 : if (eta(i,j,k) < cutoff)
530 : {
531 0 : rho_new(i,j,k,0) = rho_solid(i,j,k,0);
532 0 : M_new(i,j,k,0) = M_solid(i,j,k,0);
533 0 : M_new(i,j,k,1) = M_solid(i,j,k,1);
534 0 : E_new(i,j,k,0) = E_solid(i,j,k,0);
535 : }
536 :
537 :
538 : //Set::Vector grad_ux = Numeric::Gradient(v, i, j, k, 0, DX);
539 : //Set::Vector grad_uy = Numeric::Gradient(v, i, j, k, 1, DX);
540 :
541 6272000 : *dt_max_handle = std::fabs(cfl * DX[0] / (u(0)*eta(i,j,k) + small));
542 12544000 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl * DX[1] / (u(1)*eta(i,j,k) + small)));
543 12544000 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[0]*DX[0] / (Source(i,j,k,1)+small)));
544 12544000 : *dt_max_handle = std::min(*dt_max_handle, std::fabs(cfl_v * DX[1]*DX[1] / (Source(i,j,k,2)+small)));
545 :
546 : // Compute vorticity
547 12544000 : omega(i, j, k) = eta(i, j, k) * (gradu(1,0) - gradu(0,1));
548 :
549 6272000 : });
550 :
551 3250 : }
552 3250 : this->DynamicTimestep_SyncTimeStep(lev,dt_max);
553 :
554 3250 : }//end Advance
555 :
556 0 : void Hydro::Regrid(int lev, Set::Scalar /* time */)
557 : {
558 : BL_PROFILE("Integrator::Hydro::Regrid");
559 0 : Source_mf[lev]->setVal(0.0);
560 0 : if (lev < finest_level) return;
561 :
562 0 : Util::Message(INFO, "Regridding on level", lev);
563 : }//end regrid
564 :
565 : //void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray &a_tags, Set::Scalar time, int ngrow)
566 0 : void Hydro::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar, int)
567 : {
568 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
569 :
570 0 : const Set::Scalar* DX = geom[lev].CellSize();
571 0 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
572 :
573 : // Eta criterion for refinement
574 0 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi) {
575 0 : const amrex::Box& bx = mfi.tilebox();
576 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
577 0 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
578 :
579 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
580 0 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
581 0 : if (grad_eta.lpNorm<2>() * dr * 2 > eta_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
582 0 : });
583 0 : }
584 :
585 : // Vorticity criterion for refinement
586 0 : for (amrex::MFIter mfi(*vorticity_mf[lev], true); mfi.isValid(); ++mfi) {
587 0 : const amrex::Box& bx = mfi.tilebox();
588 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
589 0 : amrex::Array4<const Set::Scalar> const& omega = (*vorticity_mf[lev]).array(mfi);
590 :
591 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
592 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
593 0 : Set::Vector grad_omega = Numeric::Gradient(omega, i, j, k, 0, DX, sten);
594 0 : if (grad_omega.lpNorm<2>() * dr * 2 > omega_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
595 0 : });
596 0 : }
597 :
598 : // Gradu criterion for refinement
599 0 : for (amrex::MFIter mfi(*velocity_mf[lev], true); mfi.isValid(); ++mfi) {
600 0 : const amrex::Box& bx = mfi.tilebox();
601 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
602 0 : amrex::Array4<const Set::Scalar> const& v = (*velocity_mf[lev]).array(mfi);
603 :
604 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
605 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
606 0 : Set::Matrix grad_u = Numeric::Gradient(v, i, j, k, DX, sten);
607 0 : if (grad_u.lpNorm<2>() * dr * 2 > gradu_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
608 0 : });
609 0 : }
610 :
611 : // Pressure criterion for refinement
612 0 : for (amrex::MFIter mfi(*pressure_mf[lev], true); mfi.isValid(); ++mfi) {
613 0 : const amrex::Box& bx = mfi.tilebox();
614 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
615 0 : amrex::Array4<const Set::Scalar> const& p = (*pressure_mf[lev]).array(mfi);
616 :
617 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
618 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
619 0 : Set::Vector grad_p = Numeric::Gradient(p, i, j, k, 0, DX, sten);
620 0 : if (grad_p.lpNorm<2>() * dr * 2 > p_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
621 0 : });
622 0 : }
623 :
624 : // Density criterion for refinement
625 0 : for (amrex::MFIter mfi(*density_mf[lev], true); mfi.isValid(); ++mfi) {
626 0 : const amrex::Box& bx = mfi.tilebox();
627 0 : amrex::Array4<char> const& tags = a_tags.array(mfi);
628 0 : amrex::Array4<const Set::Scalar> const& rho = (*density_mf[lev]).array(mfi);
629 :
630 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
631 0 : auto sten = Numeric::GetStencil(i, j, k, bx);
632 0 : Set::Vector grad_rho = Numeric::Gradient(rho, i, j, k, 0, DX, sten);
633 0 : if (grad_rho.lpNorm<2>() * dr * 2 > rho_refinement_criterion) tags(i, j, k) = amrex::TagBox::SET;
634 0 : });
635 0 : }
636 :
637 0 : }
638 :
639 : }
640 :
641 :
642 : #endif
|
