Line data Source code
1 : #include "Flame.H"
2 : #include "IO/ParmParse.H"
3 : #include "BC/Constant.H"
4 : #include "Numeric/Stencil.H"
5 : #include "IC/Laminate.H"
6 : #include "IC/Constant.H"
7 : #include "IC/PSRead.H"
8 : #include "Numeric/Function.H"
9 : #include "IC/Expression.H"
10 : #include "IC/BMP.H"
11 : #include "IC/PNG.H"
12 : #include "Base/Mechanics.H"
13 :
14 : #include <cmath>
15 :
16 : namespace Integrator
17 : {
18 :
19 3 : Flame::Flame() : Base::Mechanics<model_type>() {}
20 :
21 3 : Flame::Flame(IO::ParmParse& pp) : Flame()
22 : {
23 9 : pp_queryclass(*this);
24 3 : }
25 :
26 : // [parser]
27 : void
28 3 : Flame::Parse(Flame& value, IO::ParmParse& pp)
29 : {
30 : BL_PROFILE("Integrator::Flame::Flame()");
31 : {
32 : // Whether to include extra fields (such as mdot, etc) in the plot output
33 18 : pp_query_default("plot_field",value.plot_field,true);
34 :
35 18 : pp_query_default("timestep", value.base_time, 1.0e-4); //Simulation timestep
36 18 : pp_query_default("pf.eps", value.pf.eps, 0.0); // Burn width thickness
37 18 : pp_query_default("pf.kappa", value.pf.kappa, 0.0); // Interface energy param
38 18 : pp_query_default("pf.gamma", value.pf.gamma, 1.0); // Scaling factor for mobility
39 18 : pp_query_default("pf.lambda", value.pf.lambda, 0.0); // Chemical potential multiplier
40 18 : pp_query_default("pf.w1", value.pf.w1, 0.0); // Unburned rest energy
41 18 : pp_query_default("pf.w12", value.pf.w12, 0.0); // Barrier energy
42 18 : pp_query_default("pf.w0", value.pf.w0, 0.0); // Burned rest energy
43 18 : pp_query_default("amr.ghost_cells", value.ghost_count, 2); // number of ghost cells in all fields
44 18 : pp_query_default("geometry.x_len", value.x_len, 0.001); // Domain x length
45 15 : pp_query_default("geometry.y_len", value.y_len, 0.001); // Domain y length
46 :
47 : //value.bc_eta = new BC::Constant(1);
48 :
49 :
50 : // Boundary conditions for phase field order params
51 6 : pp.select<BC::Constant>("pf.eta.bc", value.bc_eta, 1 );
52 9 : value.RegisterNewFab(value.eta_mf, value.bc_eta, 1, value.ghost_count, "eta", true);
53 9 : value.RegisterNewFab(value.eta_old_mf, value.bc_eta, 1, value.ghost_count, "eta_old", false);
54 :
55 : // phase field initial condition
56 9 : pp.select<IC::Laminate,IC::Constant,IC::Expression,IC::BMP,IC::PNG>("pf.eta.ic",value.ic_eta,value.geom);
57 24 : pp.forbid("eta.ic","--> you must use pf.eta.ic instead"); // delete this eventually once all tests are updated
58 : }
59 :
60 : {
61 : //IO::ParmParse pp("thermal");
62 18 : pp_query_default("thermal.on", value.thermal.on, false); // Whether to use the Thermal Transport Model
63 18 : pp_query_default("elastic.on", value.elastic.on, 0); // Whether to use Neo-hookean Elastic model
64 18 : pp_query_default("thermal.bound", value.thermal.bound, 0.0); // System Initial Temperature
65 18 : pp_query_default("elastic.traction", value.elastic.traction, 0.0); // Body force
66 15 : pp_query_default("elastic.phirefinement", value.elastic.phirefinement, 1); // Phi refinement criteria
67 :
68 :
69 :
70 3 : if (value.thermal.on) {
71 6 : pp_query_required("thermal.rho_ap", value.thermal.rho_ap); // AP Density
72 6 : pp_query_required("thermal.rho_htpb", value.thermal.rho_htpb); // HTPB Density
73 6 : pp_query_required("thermal.k_ap", value.thermal.k_ap); // AP Thermal Conductivity
74 6 : pp_query_required("thermal.k_htpb", value.thermal.k_htpb); // HTPB Thermal Conductivity
75 6 : pp_query_required("thermal.cp_ap", value.thermal.cp_ap); // AP Specific Heat
76 6 : pp_query_required("thermal.cp_htpb", value.thermal.cp_htpb); //HTPB Specific Heat
77 :
78 6 : pp_query_default("thermal.q0", value.thermal.q0, 0.0); // Baseline heat flux
79 :
80 6 : pp_query_required("thermal.m_ap", value.thermal.m_ap); // AP Pre-exponential factor for Arrhenius Law
81 6 : pp_query_required("thermal.m_htpb", value.thermal.m_htpb); // HTPB Pre-exponential factor for Arrhenius Law
82 6 : pp_query_required("thermal.E_ap", value.thermal.E_ap); // AP Activation Energy for Arrhenius Law
83 6 : pp_query_required("thermal.E_htpb", value.thermal.E_htpb); // HTPB Activation Energy for Arrhenius Law
84 :
85 6 : pp_query_default("thermal.hc", value.thermal.hc, 1.0); // Used to change heat flux units
86 6 : pp_query_default("thermal.massfraction", value.thermal.massfraction, 0.8); // Systen AP mass fraction
87 6 : pp_query_default("thermal.mlocal_ap", value.thermal.mlocal_ap, 0.0); // AP mass flux reference value
88 7 : pp_query_default("thermal.mlocal_htpb", value.thermal.mlocal_htpb, 0.0); // HTPB mass flux reference value
89 6 : pp_query_default("thermal.mlocal_comb", value.thermal.mlocal_comb, 0.0); // AP/HTPB mass flux reference value
90 :
91 5 : pp_query_default("thermal.T_fluid", value.thermal.T_fluid, 300.0); // Temperature of the Standin Fluid
92 :
93 1 : pp_query("thermal.disperssion1", value.thermal.disperssion1); // K; dispersion variables are use to set the outter field properties for the void grain case.
94 1 : pp_query("thermal.disperssion2", value.thermal.disperssion2); // rho; dispersion variables are use to set the outter field properties for the void grain case.
95 1 : pp_query("thermal.disperssion3", value.thermal.disperssion3); // cp; dispersion variables are use to set the outter field properties for the void grain case.
96 :
97 6 : pp_query_default("thermal.modeling_ap", value.thermal.modeling_ap, 1.0); // Scaling factor for AP thermal conductivity (default = 1.0)
98 5 : pp_query_default("thermal.modeling_htpb", value.thermal.modeling_htpb, 1.0); // Scaling factor for HTPB thermal conductivity (default = 1.0)
99 :
100 :
101 : //Temperature boundary condition
102 2 : pp.select_default<BC::Constant>("thermal.temp.bc", value.bc_temp, 1);
103 :
104 3 : value.RegisterNewFab(value.temp_mf, value.bc_temp, 1, value.ghost_count + 1, "temp", true);
105 3 : value.RegisterNewFab(value.temp_old_mf, value.bc_temp, 1, value.ghost_count + 1, "temp_old", false);
106 3 : value.RegisterNewFab(value.temps_mf, value.bc_temp, 1, value.ghost_count + 1, "temps", false);
107 3 : value.RegisterNewFab(value.temps_old_mf, value.bc_temp, 1, value.ghost_count + 1, "temps_old", false);
108 :
109 3 : value.RegisterNewFab(value.mdot_mf, value.bc_eta, 1, value.ghost_count + 1, "mdot", value.plot_field);
110 3 : value.RegisterNewFab(value.mob_mf, value.bc_eta, 1, value.ghost_count + 1, "mob", value.plot_field);
111 3 : value.RegisterNewFab(value.alpha_mf, value.bc_temp, 1, value.ghost_count + 1, "alpha", value.plot_field);
112 3 : value.RegisterNewFab(value.heatflux_mf, value.bc_temp, 1, value.ghost_count + 1, "heatflux", value.plot_field);
113 3 : value.RegisterNewFab(value.laser_mf, value.bc_temp, 1, value.ghost_count + 1, "laser", value.plot_field);
114 :
115 2 : value.RegisterIntegratedVariable(&value.volume, "total_area");
116 2 : value.RegisterIntegratedVariable(&value.area, "Interface_area");
117 2 : value.RegisterIntegratedVariable(&value.chamber_area, "chamber_area", false);
118 2 : value.RegisterIntegratedVariable(&value.massflux, "mass_flux");
119 2 : value.RegisterIntegratedVariable(&value.chamber_pressure, "Pressure", false);
120 :
121 : // laser initial condition
122 2 : pp.select_default<IC::Constant,IC::Expression>("laser.ic",value.ic_laser, value.geom);
123 :
124 : // thermal initial condition
125 3 : pp.select_default<IC::Constant,IC::Expression,IC::BMP,IC::PNG>("temp.ic",value.thermal.ic_temp,value.geom);
126 : }
127 : }
128 :
129 : {
130 15 : pp_query_default("pressure.P", value.pressure.P, 1.0); // Constant pressure value
131 3 : if (value.thermal.on)
132 : {
133 6 : pp_query_required("pressure.a1", value.pressure.arrhenius.a1); // Surgate heat flux model paramater - AP
134 6 : pp_query_required("pressure.a2", value.pressure.arrhenius.a2); // Surgate heat flux model paramater - HTPB
135 6 : pp_query_required("pressure.a3", value.pressure.arrhenius.a3); // Surgate heat flux model paramater - Total
136 6 : pp_query_required("pressure.b1", value.pressure.arrhenius.b1); // Surgate heat flux model paramater - AP
137 6 : pp_query_required("pressure.b2", value.pressure.arrhenius.b2); // Surgate heat flux model paramater - HTPB
138 6 : pp_query_required("pressure.b3", value.pressure.arrhenius.b3); // Surgate heat flux model paramater - Total
139 6 : pp_query_required("pressure.c1", value.pressure.arrhenius.c1); // Surgate heat flux model paramater - Total
140 6 : pp_query_default("pressure.mob_ap", value.pressure.arrhenius.mob_ap, 0); // Whether to include pressure to the arrhenius law
141 6 : pp_query_default("pressure.dependency", value.pressure.arrhenius.dependency, 1); // Whether to use pressure to determined the reference Zeta
142 6 : pp_query_default("pressure.h1", value.pressure.arrhenius.h1, 1.81); // Surgate heat flux model paramater - Homogenized
143 6 : pp_query_default("pressure.h2", value.pressure.arrhenius.h2, 1.34); // Surgate heat flux model paramater - Homogenized
144 :
145 : }
146 : else
147 : {
148 12 : pp_query_required("pressure.r_ap", value.pressure.power.r_ap); // AP power pressure law parameter (r*P^n)
149 12 : pp_query_required("pressure.r_htpb", value.pressure.power.r_htpb); // HTPB power pressure law parameter (r*P^n)
150 12 : pp_query_required("pressure.r_comb", value.pressure.power.r_comb); // AP/HTPB power pressure law parameter (r*P^n)
151 12 : pp_query_required("pressure.n_ap", value.pressure.power.n_ap); // AP power pressure law parameter (r*P^n)
152 12 : pp_query_required("pressure.n_htpb", value.pressure.power.n_htpb); // HTPB power pressure law parameter (r*P^n)
153 12 : pp_query_required("pressure.n_comb", value.pressure.power.n_comb); // AP/HTPB power pressure law parameter (r*P^n)
154 : }
155 18 : pp_query_default("variable_pressure", value.variable_pressure, 0); // Whether to compute the pressure evolution
156 18 : pp_query_default("homogeneousSystem", value.homogeneousSystem, 0); // Whether to initialize Phi with homogenized properties
157 : }
158 :
159 18 : pp_query_default("amr.refinement_criterion", value.m_refinement_criterion, 0.001);// Refinement criterion for eta field
160 18 : pp_query_default("amr.refinement_criterion_temp", value.t_refinement_criterion, 0.001);// Refinement criterion for temperature field
161 18 : pp_query_default("amr.refinament_restriction", value.t_refinement_restriction, 0.1);// Eta value to restrict the refinament for the temperature field
162 18 : pp_query_default("amr.phi_refinement_criterion", value.phi_refinement_criterion, 1.0e100);// Refinement criterion for phi field [infinity]
163 18 : pp_query_default("small", value.small, 1.0e-8); // Lowest value of Eta.
164 :
165 : {
166 : // The material field is referred to as :math:`\phi(\mathbf{x})` and is
167 : // specified using these parameters.
168 6 : std::string phi_ic_type = "packedspheres";
169 21 : pp_query_validate("phi.ic.type", phi_ic_type, {"psread", "laminate", "expression", "constant", "bmp", "png"}); // IC type (psread, laminate, constant)
170 3 : if (phi_ic_type == "psread") {
171 0 : value.ic_phi = new IC::PSRead(value.geom, pp, "phi.ic.psread");
172 : //value.ic_phicell = new IC::PSRead(value.geom, pp, "phi.ic.psread");
173 0 : pp_query_default("phi.ic.psread.eps", value.zeta, 1.0e-5); // AP/HTPB interface length
174 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
175 : }
176 3 : else if (phi_ic_type == "laminate") {
177 9 : value.ic_phi = new IC::Laminate(value.geom, pp, "phi.ic.laminate");
178 : //value.ic_phicell = new IC::Laminate(value.geom, pp, "phi.ic.laminate");
179 18 : pp_query_default("phi.ic.laminate.eps", value.zeta, 1.0e-5); // AP/HTPB interface length
180 18 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
181 : }
182 0 : else if (phi_ic_type == "expression") {
183 0 : value.ic_phi = new IC::Expression(value.geom, pp, "phi.ic.expression");
184 : //value.ic_phicell = new IC::Expression(value.geom, pp, "phi.ic.expression");
185 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
186 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
187 : }
188 0 : else if (phi_ic_type == "constant") {
189 0 : value.ic_phi = new IC::Constant(value.geom, pp, "phi.ic.constant");
190 : //value.ic_phicell = new IC::Constant(value.geom, pp, "phi.ic.constant");
191 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
192 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
193 : }
194 0 : else if (phi_ic_type == "bmp") {
195 0 : value.ic_phi = new IC::BMP(value.geom, pp, "phi.ic.bmp");
196 : //value.ic_phicell = new IC::BMP(value.geom, pp, "phi.ic.bmp");
197 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
198 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
199 : }
200 0 : else if (phi_ic_type == "png") {
201 0 : value.ic_phi = new IC::PNG(value.geom, pp, "phi.ic.png");
202 0 : pp_query_default("phi.zeta_0", value.zeta_0, 1.0e-5); // Reference interface length for heat integration
203 0 : pp_query_default("phi.zeta", value.zeta, 1.0e-5); // AP/HTPB interface length
204 : }
205 0 : else Util::Abort(INFO, "Invalid IC type ", phi_ic_type);
206 : //value.RegisterNewFab(value.phi_mf, 1, "phi_cell", true);
207 9 : value.RegisterNodalFab(value.phi_mf, 1, value.ghost_count + 1, "phi", true);
208 : //value.RegisterNewFab(value.phicell_mf, value.bc_eta, 1, value.ghost_count + 1, "phi", true);
209 3 : }
210 :
211 15 : pp.queryclass<Base::Mechanics<model_type>>("elastic",value);
212 :
213 3 : if (value.m_type != Type::Disable)
214 : {
215 0 : value.elastic.Tref = value.thermal.bound;
216 0 : pp_query_default("Tref", value.elastic.Tref, 300.0); // Initial temperature for thermal expansion computation
217 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_ap", value.elastic.model_ap);
218 0 : pp.queryclass<Model::Solid::Finite::NeoHookeanPredeformed>("model_htpb", value.elastic.model_htpb);
219 :
220 0 : value.bc_psi = new BC::Nothing();
221 0 : value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, value.ghost_count, "psi", value.plot_psi);
222 0 : value.psi_on = true;
223 : }
224 :
225 : bool allow_unused;
226 : // Set this to true to allow unused inputs without error.
227 : // (Not recommended.)
228 15 : pp.query_default("allow_unused",allow_unused,false);
229 3 : if (!allow_unused && pp.AnyUnusedInputs())
230 : {
231 0 : Util::Warning(INFO,"The following inputs were specified but not used:");
232 0 : pp.AllUnusedInputs();
233 0 : Util::Exception(INFO,"Aborting. Specify 'allow_unused=True` to ignore this error.");
234 : }
235 :
236 3 : }
237 :
238 24 : void Flame::Initialize(int lev)
239 : {
240 : BL_PROFILE("Integrator::Flame::Initialize");
241 24 : Base::Mechanics<model_type>::Initialize(lev);
242 :
243 24 : ic_eta->Initialize(lev, eta_mf);
244 24 : ic_eta->Initialize(lev, eta_old_mf);
245 24 : ic_phi->Initialize(lev, phi_mf);
246 : //ic_phicell->Initialize(lev, phicell_mf);
247 :
248 24 : if (elastic.on) {
249 0 : psi_mf[lev]->setVal(1.0);
250 0 : rhs_mf[lev]->setVal(Set::Vector::Zero());
251 : }
252 24 : if (thermal.on) {
253 8 : if (thermal.ic_temp)
254 : {
255 8 : thermal.ic_temp->Initialize(lev,temp_mf);
256 8 : thermal.ic_temp->Initialize(lev,temp_old_mf);
257 8 : thermal.ic_temp->Initialize(lev,temps_mf);
258 8 : thermal.ic_temp->Initialize(lev,temps_old_mf);
259 : }
260 : else
261 : {
262 0 : temp_mf[lev]->setVal(thermal.bound);
263 0 : temp_old_mf[lev]->setVal(thermal.bound);
264 0 : temps_mf[lev]->setVal(thermal.bound);
265 0 : temps_old_mf[lev]->setVal(thermal.bound);
266 : }
267 8 : alpha_mf[lev]->setVal(0.0);
268 8 : mob_mf[lev]->setVal(0.0);
269 8 : mdot_mf[lev]->setVal(0.0);
270 8 : heatflux_mf[lev]->setVal(0.0);
271 : // pressure_mf[lev]->setVal(1.0); // [error]
272 8 : thermal.w1 = 0.2 * pressure.P + 0.9;
273 8 : thermal.T_fluid = thermal.bound;
274 8 : ic_laser->Initialize(lev, laser_mf);
275 8 : thermal.mlocal_htpb = 685000.0 - 850e3 * thermal.massfraction;
276 : }
277 24 : if (variable_pressure) pressure.P = 1.0;
278 24 : }
279 :
280 0 : void Flame::UpdateModel(int /*a_step*/, Set::Scalar /*a_time*/)
281 : {
282 0 : if (m_type == Base::Mechanics<model_type>::Type::Disable) return;
283 :
284 0 : for (int lev = 0; lev <= finest_level; ++lev)
285 : {
286 0 : amrex::Box domain = this->geom[lev].Domain();
287 0 : domain.convert(amrex::IntVect::TheNodeVector());
288 0 : const Set::Scalar* DX = geom[lev].CellSize();
289 :
290 : //psi_mf[lev]->setVal(1.0);
291 0 : phi_mf[lev]->FillBoundary();
292 : //phicell_mf[lev]->FillBoundary();
293 0 : eta_mf[lev]->FillBoundary();
294 0 : temp_mf[lev]->FillBoundary();
295 :
296 0 : for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
297 : {
298 0 : amrex::Box smallbox = mfi.nodaltilebox();
299 0 : amrex::Box bx = mfi.grownnodaltilebox() & domain;
300 : //amrex::Box bx = mfi.nodaltilebox();
301 : //bx.grow(1);
302 0 : amrex::Array4<model_type> const& model = model_mf[lev]->array(mfi);
303 0 : amrex::Array4<const Set::Scalar> const& phi = phi_mf[lev]->array(mfi);
304 0 : amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
305 0 : amrex::Array4<Set::Vector> const& rhs = rhs_mf[lev]->array(mfi);
306 : // amrex::Array4<const Set::Scalar> const& Pressure = pressure_mf[lev]->array(mfi); // [error]
307 :
308 0 : if (elastic.on)
309 : {
310 0 : amrex::Array4<const Set::Scalar> const& temp = temp_mf[lev]->array(mfi);
311 0 : amrex::ParallelFor(smallbox, [=] AMREX_GPU_DEVICE(int i, int j, int k)
312 : {
313 0 : Set::Vector grad_eta = Numeric::CellGradientOnNode(eta, i, j, k, 0, DX);
314 0 : rhs(i, j, k) = elastic.traction * grad_eta;
315 0 : });
316 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
317 : {
318 0 : Set::Scalar phi_avg = phi(i, j, k, 0);
319 0 : Set::Scalar temp_avg = Numeric::Interpolate::CellToNodeAverage(temp, i, j, k, 0);
320 0 : model_type model_ap = elastic.model_ap;
321 0 : model_ap.F0 -= Set::Matrix::Identity();
322 0 : model_ap.F0 *= (temp_avg - elastic.Tref);
323 0 : model_ap.F0 += Set::Matrix::Identity();
324 0 : model_type model_htpb = elastic.model_htpb;
325 0 : model_htpb.F0 -= Set::Matrix::Identity();
326 0 : model_htpb.F0 *= (temp_avg - elastic.Tref);
327 0 : model_htpb.F0 += Set::Matrix::Identity();
328 :
329 0 : model(i, j, k) = model_ap * phi_avg + model_htpb * (1. - phi_avg);
330 0 : });
331 : }
332 : else
333 : {
334 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
335 : {
336 0 : Set::Scalar phi_avg = Numeric::Interpolate::CellToNodeAverage(phi, i, j, k, 0);
337 : //phi_avg = phi(i,j,k,0);
338 0 : model_type model_ap = elastic.model_ap;
339 0 : model_ap.F0 *= Set::Matrix::Zero();
340 0 : model_type model_htpb = elastic.model_htpb;
341 0 : model_htpb.F0 *= Set::Matrix::Zero();
342 0 : model(i, j, k) = model_ap * phi_avg + model_htpb * (1. - phi_avg);
343 0 : });
344 : }
345 0 : }
346 0 : Util::RealFillBoundary(*model_mf[lev], geom[lev]);
347 :
348 0 : psi_mf[lev]->setVal(1.0);
349 0 : amrex::MultiFab::Copy(*psi_mf[lev], *eta_mf[lev], 0, 0, 1, psi_mf[lev]->nGrow());
350 : }
351 : }
352 :
353 212 : void Flame::TimeStepBegin(Set::Scalar a_time, int a_iter)
354 : {
355 : BL_PROFILE("Integrator::Flame::TimeStepBegin");
356 212 : Base::Mechanics<model_type>::TimeStepBegin(a_time, a_iter);
357 212 : if (thermal.on) {
358 90 : for (int lev = 0; lev <= finest_level; ++lev)
359 80 : ic_laser->Initialize(lev, laser_mf, a_time);
360 : }
361 212 : }
362 :
363 212 : void Flame::TimeStepComplete(Set::Scalar /*a_time*/, int /*a_iter*/)
364 : {
365 : BL_PROFILE("Integrator::Flame::TimeStepComplete");
366 212 : if (variable_pressure) {
367 0 : Set::Scalar domain_area = x_len * y_len;
368 0 : chamber_pressure = pressure.P;
369 0 : chamber_area = domain_area - volume;
370 0 : Util::Message(INFO, "Mass = ", massflux);
371 0 : Util::Message(INFO, "Pressure = ", pressure.P);
372 : }
373 212 : }
374 :
375 54060 : void Flame::Advance(int lev, Set::Scalar time, Set::Scalar dt)
376 : {
377 : BL_PROFILE("Integrador::Flame::Advance");
378 54060 : Base::Mechanics<model_type>::Advance(lev, time, dt);
379 54060 : const Set::Scalar* DX = geom[lev].CellSize();
380 :
381 : if (true) //lev == finest_level) //(true)
382 : {
383 54060 : if (thermal.on)
384 : {
385 2550 : std::swap(eta_old_mf[lev], eta_mf[lev]);
386 2550 : std::swap(temp_old_mf[lev], temp_mf[lev]);
387 2550 : std::swap(temps_old_mf[lev], temps_mf[lev]);
388 :
389 : Numeric::Function::Polynomial<4> w(pf.w0,
390 : 0.0,
391 2550 : -5.0 * pf.w1 + 16.0 * pf.w12 - 11.0 * pf.w0,
392 2550 : 14.0 * pf.w1 - 32.0 * pf.w12 + 18.0 * pf.w0,
393 2550 : -8.0 * pf.w1 + 16.0 * pf.w12 - 8.0 * pf.w0);
394 2550 : Numeric::Function::Polynomial<3> dw = w.D();
395 :
396 13100 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
397 : {
398 10550 : const amrex::Box& bx = mfi.tilebox();
399 : // Phase fields
400 10550 : amrex::Array4<Set::Scalar> const& etanew = (*eta_mf[lev]).array(mfi);
401 10550 : amrex::Array4<const Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
402 10550 : amrex::Array4<const Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
403 : //amrex::Array4<const Set::Scalar> const& phicell = (*phicell_mf[lev]).array(mfi);
404 : // Heat transfer fields
405 10550 : amrex::Array4<Set::Scalar> const& tempnew = (*temp_mf[lev]).array(mfi);
406 10550 : amrex::Array4<const Set::Scalar> const& temp = (*temp_old_mf[lev]).array(mfi);
407 10550 : amrex::Array4<Set::Scalar> const& alpha = (*alpha_mf[lev]).array(mfi);
408 10550 : amrex::Array4<Set::Scalar> const& tempsnew = (*temps_mf[lev]).array(mfi);
409 10550 : amrex::Array4<Set::Scalar> const& temps = (*temps_old_mf[lev]).array(mfi);
410 10550 : amrex::Array4<Set::Scalar> const& laser = (*laser_mf[lev]).array(mfi);
411 : // Diagnostic fields
412 10550 : amrex::Array4<Set::Scalar> const& mob = (*mob_mf[lev]).array(mfi);
413 10550 : amrex::Array4<Set::Scalar> const& mdot = (*mdot_mf[lev]).array(mfi);
414 10550 : amrex::Array4<Set::Scalar> const& heatflux = (*heatflux_mf[lev]).array(mfi);
415 :
416 10550 : if (variable_pressure) {
417 0 : pressure.P = exp(0.00075 * massflux);
418 0 : if (pressure.P > 10.0) {
419 0 : pressure.P = 10.0;
420 : }
421 0 : else if (pressure.P <= 0.99) {
422 0 : pressure.P = 0.99;
423 : }
424 0 : elastic.traction = pressure.P;
425 :
426 : }
427 :
428 10550 : Set::Scalar zeta_2 = 0.000045 - pressure.P * 6.42e-6;
429 : Set::Scalar zeta_1;
430 10550 : if (pressure.arrhenius.dependency == 1) zeta_1 = zeta_2;
431 0 : else zeta_1 = zeta_0;
432 :
433 10550 : Set::Scalar k1 = pressure.arrhenius.a1 * pressure.P + pressure.arrhenius.b1 - zeta_1 / zeta;
434 10550 : Set::Scalar k2 = pressure.arrhenius.a2 * pressure.P + pressure.arrhenius.b2 - zeta_1 / zeta;
435 10550 : Set::Scalar k3 = 4.0 * log((pressure.arrhenius.c1 * pressure.P * pressure.P + pressure.arrhenius.a3 * pressure.P + pressure.arrhenius.b3) - k1 / 2.0 - k2 / 2.0);
436 10550 : Set::Scalar k4 = pressure.arrhenius.h1 * pressure.P + pressure.arrhenius.h2;
437 :
438 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
439 : {
440 1747840 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
441 1747840 : Set::Scalar eta_lap = Numeric::Laplacian(eta, i, j, k, 0, DX);
442 : Set::Scalar K; // Calculate effective thermal conductivity
443 : Set::Scalar rho; // No special interface mixure rule is needed here.
444 : Set::Scalar cp;
445 1747840 : if (homogeneousSystem) {
446 0 : K = (thermal.modeling_ap * thermal.k_ap * thermal.massfraction + thermal.modeling_htpb * thermal.k_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion1 * (1. - phi_avg); // Calculate effective thermal conductivity
447 0 : rho = (thermal.rho_ap * thermal.massfraction + thermal.rho_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion2 * (1. - phi_avg); // No special interface mixure rule is needed here.
448 0 : cp = (thermal.cp_ap * thermal.massfraction + thermal.cp_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion3 * (1. - phi_avg);
449 : }
450 : else {
451 1747840 : K = thermal.modeling_ap * thermal.k_ap * phi_avg + thermal.modeling_htpb * thermal.k_htpb * (1.0 - phi_avg); // Calculate effective thermal conductivity
452 1747840 : rho = thermal.rho_ap * phi_avg + thermal.rho_htpb * (1.0 - phi_avg); // No special interface mixure rule is needed here.
453 1747840 : cp = thermal.cp_ap * phi_avg + thermal.cp_htpb * (1.0 - phi_avg);
454 : }
455 3495680 : Set::Scalar df_deta = ((pf.lambda / pf.eps) * dw(eta(i, j, k)) - pf.eps * pf.kappa * eta_lap);
456 5243520 : etanew(i, j, k) = eta(i, j, k) - mob(i, j, k) * dt * df_deta;
457 3495680 : if (etanew(i, j, k) <= small) etanew(i, j, k) = small;
458 :
459 1747840 : alpha(i, j, k) = K / rho / cp; // Calculate thermal diffusivity and store in fiel
460 5243520 : mdot(i, j, k) = rho * fabs(eta(i, j, k) - etanew(i, j, k)) / dt; // deta/dt
461 :
462 6991360 : if (isnan(etanew(i, j, k)) || isnan(alpha(i, j, k)) || isnan(mdot(i, j, k))) {
463 0 : Util::Message(INFO, etanew(i, j, k), "etanew contains nan (i=", i, " j= ", j, ")");
464 0 : Util::Message(INFO, mdot(i, j, k), "mdot contains nan (i=", i, " j= ", j, ")");
465 0 : Util::Message(INFO, alpha(i, j, k), "alpha contains nan (i=", i, " j= ", j, ")");
466 0 : Util::Abort(INFO);
467 : }
468 1747840 : });
469 :
470 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
471 : {
472 1747840 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
473 : //phi_avg = phicell(i,j,k);
474 : Set::Scalar K;
475 : Set::Scalar qflux;
476 : Set::Scalar mlocal;
477 : Set::Scalar mdota;
478 : Set::Scalar mbase;
479 :
480 1747840 : if (homogeneousSystem) {
481 0 : qflux = k4 * phi_avg;
482 0 : mlocal = (thermal.mlocal_ap) * thermal.massfraction + (thermal.mlocal_htpb) * (1.0 - thermal.massfraction);
483 0 : mdota = fabs(mdot(i, j, k));
484 0 : mbase = tanh(4.0 * mdota / (mlocal));
485 0 : K = (thermal.modeling_ap * thermal.k_ap * thermal.massfraction + thermal.modeling_htpb * thermal.k_htpb * (1.0 - thermal.massfraction)) * phi_avg + thermal.disperssion1 * (1. - phi_avg);
486 0 : heatflux(i, j, k) = (laser(i, j, k) * phi_avg + thermal.hc * mbase * qflux) / K;
487 : }
488 : else {
489 1747840 : qflux = k1 * phi_avg +
490 1747840 : k2 * (1.0 - phi_avg) +
491 1747840 : (zeta_1 / zeta) * exp(k3 * phi_avg * (1.0 - phi_avg));
492 1747840 : mlocal = (thermal.mlocal_ap) * phi_avg + (thermal.mlocal_htpb) * (1.0 - phi_avg) + thermal.mlocal_comb * phi_avg * (1.0 - phi_avg);
493 1747840 : mdota = fabs(mdot(i, j, k));
494 1747840 : mbase = tanh(4.0 * mdota / (mlocal));
495 :
496 1747840 : K = thermal.modeling_ap * thermal.k_ap * phi_avg + thermal.modeling_htpb * thermal.k_htpb * (1.0 - phi_avg);
497 5243520 : heatflux(i, j, k) = (thermal.hc * mbase * qflux + laser(i, j, k)) / K;
498 : }
499 :
500 3495680 : if (isnan(heatflux(i, j, k))) {
501 0 : Util::Message(INFO, heatflux(i, j, k), "heat contains nan (i=", i, " j= ", j, ")");
502 0 : Util::Abort(INFO);
503 : }
504 1747840 : });
505 :
506 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
507 : {
508 1747840 : auto sten = Numeric::GetStencil(i, j, k, bx);
509 1747840 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
510 1747840 : Set::Vector grad_temp = Numeric::Gradient(temp, i, j, k, 0, DX);
511 1747840 : Set::Scalar lap_temp = Numeric::Laplacian(temp, i, j, k, 0, DX);
512 1747840 : Set::Scalar grad_eta_mag = grad_eta.lpNorm<2>();
513 1747840 : Set::Vector grad_alpha = Numeric::Gradient(alpha, i, j, k, 0, DX, sten);
514 1747840 : Set::Scalar dTdt = 0.0;
515 3495680 : dTdt += grad_eta.dot(grad_temp * alpha(i, j, k));
516 3495680 : dTdt += grad_alpha.dot(eta(i, j, k) * grad_temp);
517 3495680 : dTdt += eta(i, j, k) * alpha(i, j, k) * lap_temp;
518 3495680 : dTdt += alpha(i, j, k) * heatflux(i, j, k) * grad_eta_mag;
519 : Set::Scalar Tsolid;
520 5243520 : Tsolid = dTdt + temps(i, j, k) * (etanew(i, j, k) - eta(i, j, k)) / dt;
521 3495680 : tempsnew(i, j, k) = temps(i, j, k) + dt * Tsolid;
522 6991360 : tempnew(i, j, k) = etanew(i, j, k) * tempsnew(i, j, k) + (1.0 - etanew(i, j, k)) * thermal.T_fluid;
523 5243520 : if (isnan(tempsnew(i, j, k)) || isnan(temps(i, j, k))) {
524 0 : Util::Message(INFO, tempsnew(i, j, k), "tempsnew contains nan (i=", i, " j= ", j, ")");
525 0 : Util::Message(INFO, temps(i, j, k), "temps contains nan (i=", i, " j= ", j, ")");
526 0 : Util::Abort(INFO);
527 : }
528 :
529 1747840 : });
530 :
531 10550 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
532 : {
533 1747840 : Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
534 : //phi_avg = phicell(i,j,k);
535 :
536 : Set::Scalar L;
537 1747840 : if (pressure.arrhenius.mob_ap == 1) L = thermal.m_ap * pressure.P * exp(-thermal.E_ap / tempnew(i, j, k)) * phi_avg;
538 3495680 : else L = thermal.m_ap * exp(-thermal.E_ap / tempnew(i, j, k)) * phi_avg;
539 1747840 : L += thermal.m_htpb * exp(-thermal.E_htpb / tempnew(i, j, k)) * (1.0 - phi_avg);
540 : //L += thermal.m_comb * (0.5 * tempnew(i, j, k) / thermal.bound) * phi(i, j, k) * (1.0 - phi(i, j, k));
541 3501038 : if (tempnew(i, j, k) <= thermal.bound) mob(i, j, k) = 0;
542 3484964 : else mob(i, j, k) = L;
543 : //mob(i,j,k) = L;
544 3495680 : if (isnan(mob(i, j, k))) {
545 0 : Util::Message(INFO, mob(i, j, k), "mob contains nan (i=", i, " j= ", j, ")");
546 0 : Util::Abort(INFO);
547 : }
548 1747840 : });
549 2550 : } // MFi For loop
550 :
551 : } // thermal IF
552 : else
553 : {
554 51510 : std::swap(eta_old_mf[lev], eta_mf[lev]);
555 : Numeric::Function::Polynomial<4> w(pf.w0,
556 : 0.0,
557 51510 : -5.0 * pf.w1 + 16.0 * pf.w12 - 11.0 * pf.w0,
558 51510 : 14.0 * pf.w1 - 32.0 * pf.w12 + 18.0 * pf.w0,
559 51510 : -8.0 * pf.w1 + 16.0 * pf.w12 - 8.0 * pf.w0
560 51510 : );
561 51510 : Numeric::Function::Polynomial<3> dw = w.D();
562 :
563 149676 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
564 : {
565 98166 : const amrex::Box& bx = mfi.tilebox();
566 : //Phase Fields
567 98166 : amrex::Array4<Set::Scalar> const& etanew = (*eta_mf[lev]).array(mfi);
568 98166 : amrex::Array4<Set::Scalar> const& eta = (*eta_old_mf[lev]).array(mfi);
569 98166 : amrex::Array4<Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
570 :
571 98166 : Set::Scalar fmod_ap = pressure.power.r_ap * pow(pressure.P, pressure.power.n_ap);
572 98166 : Set::Scalar fmod_htpb = pressure.power.r_htpb * pow(pressure.P, pressure.power.n_htpb);
573 98166 : Set::Scalar fmod_comb = pressure.power.r_comb * pow(pressure.P, pressure.power.n_comb);
574 :
575 98166 : pressure.power.a_fit = -1.16582 * sin(pressure.P) - 0.681788 * cos(pressure.P) + 3.3563;
576 98166 : pressure.power.b_fit = -0.708225 * sin(pressure.P) + 0.548067 * cos(pressure.P) + 1.55985;
577 98166 : pressure.power.c_fit = -0.0130849 * sin(pressure.P) - 0.03597 * cos(pressure.P) + 0.00725694;
578 :
579 :
580 98166 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
581 : {
582 : //Set::Scalar phi_avg = Numeric::Interpolate::NodeToCellAverage(phi, i, j, k, 0);
583 : Set::Scalar L;
584 95048768 : if (homogeneousSystem == 1) {
585 0 : Set::Vector grad_eta = Numeric::Gradient(eta, i, j, k, 0, DX);
586 0 : Set::Scalar angle = acos(grad_eta[0] / grad_eta.lpNorm<2>()) * 180 / 3.1415;
587 0 : if (angle > 90) angle = angle - 90.0;
588 0 : if (angle > 180) angle = angle - 180.0;
589 0 : if (angle > 270) angle = angle - 270.0;
590 0 : L = pressure.power.a_fit + pressure.power.b_fit * exp(-pressure.power.c_fit * angle);
591 : }
592 : else {
593 : Set::Scalar fs_actual;
594 95048768 : fs_actual = fmod_ap * phi(i, j, k)
595 95048768 : + fmod_htpb * (1.0 - phi(i, j, k))
596 190097536 : + 4.0 * fmod_comb * phi(i, j, k) * (1.0 - phi(i, j, k));
597 95048768 : L = fs_actual / pf.gamma / (pf.w1 - pf.w0);
598 : }
599 95048768 : Set::Scalar eta_lap = Numeric::Laplacian(eta, i, j, k, 0, DX);
600 190097536 : Set::Scalar df_deta = ((pf.lambda / pf.eps) * dw(eta(i, j, k)) - pf.eps * pf.kappa * eta_lap);
601 190097536 : etanew(i, j, k) = eta(i, j, k) - L * dt * df_deta;
602 :
603 308123190 : if (etanew(i, j, k) > eta(i, j, k)) etanew(i, j, k) = eta(i, j, k);
604 95048768 : });
605 51510 : } // MFI For Loop
606 : } // thermal ELSE
607 : }// First IF
608 54060 : } //Function
609 :
610 :
611 263 : void Flame::TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar time, int ngrow)
612 : {
613 : BL_PROFILE("Integrator::Flame::TagCellsForRefinement");
614 263 : Base::Mechanics<model_type>::TagCellsForRefinement(lev, a_tags, time, ngrow);
615 :
616 263 : const Set::Scalar* DX = geom[lev].CellSize();
617 263 : Set::Scalar dr = sqrt(AMREX_D_TERM(DX[0] * DX[0], +DX[1] * DX[1], +DX[2] * DX[2]));
618 :
619 : // Eta criterion for refinement
620 709 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
621 : {
622 446 : const amrex::Box& bx = mfi.tilebox();
623 446 : amrex::Array4<char> const& tags = a_tags.array(mfi);
624 446 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
625 :
626 446 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
627 : {
628 78848 : Set::Vector gradeta = Numeric::Gradient(eta, i, j, k, 0, DX);
629 103225 : if (gradeta.lpNorm<2>() * dr * 2 > m_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
630 47698 : tags(i, j, k) = amrex::TagBox::SET;
631 78848 : });
632 263 : }
633 :
634 : // Phi criterion for refinement
635 263 : if (elastic.phirefinement) {
636 709 : for (amrex::MFIter mfi(*eta_mf[lev], true); mfi.isValid(); ++mfi)
637 : {
638 446 : const amrex::Box& bx = mfi.tilebox();
639 446 : amrex::Array4<char> const& tags = a_tags.array(mfi);
640 446 : amrex::Array4<const Set::Scalar> const& phi = (*phi_mf[lev]).array(mfi);
641 :
642 446 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
643 : {
644 78848 : Set::Vector gradphi = Numeric::Gradient(phi, i, j, k, 0, DX);
645 78848 : if (gradphi.lpNorm<2>() * dr >= phi_refinement_criterion)
646 0 : tags(i, j, k) = amrex::TagBox::SET;
647 78848 : });
648 263 : }
649 : }
650 :
651 :
652 : // Thermal criterion for refinement
653 263 : if (thermal.on) {
654 507 : for (amrex::MFIter mfi(*temp_mf[lev], true); mfi.isValid(); ++mfi)
655 : {
656 345 : const amrex::Box& bx = mfi.tilebox();
657 345 : amrex::Array4<char> const& tags = a_tags.array(mfi);
658 345 : amrex::Array4<const Set::Scalar> const& temp = (*temp_mf[lev]).array(mfi);
659 345 : amrex::Array4<const Set::Scalar> const& eta = (*eta_mf[lev]).array(mfi);
660 345 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
661 : {
662 50928 : Set::Vector tempgrad = Numeric::Gradient(temp, i, j, k, 0, DX);
663 50928 : if (tempgrad.lpNorm<2>() * dr > t_refinement_criterion && eta(i, j, k) >= t_refinement_restriction)
664 0 : tags(i, j, k) = amrex::TagBox::SET;
665 50928 : });
666 162 : }
667 : }
668 :
669 :
670 263 : }
671 :
672 15 : void Flame::Regrid(int lev, Set::Scalar time)
673 : {
674 : BL_PROFILE("Integrator::Flame::Regrid");
675 : //if (lev < finest_level) return;
676 : //phi_mf[lev]->setVal(0.0);
677 15 : ic_phi->Initialize(lev, phi_mf, time);
678 : //ic_phicell->Initialize(lev, phi_mf, time);
679 15 : }
680 :
681 220 : void Flame::Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
682 : const amrex::MFIter& mfi, const amrex::Box& box)
683 : {
684 : BL_PROFILE("Flame::Integrate");
685 220 : const Set::Scalar* DX = geom[amrlev].CellSize();
686 220 : Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
687 220 : amrex::Array4<amrex::Real> const& eta = (*eta_mf[amrlev]).array(mfi);
688 220 : amrex::Array4<amrex::Real> const& mdot = (*mdot_mf[amrlev]).array(mfi);
689 220 : if (variable_pressure) {
690 0 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
691 : {
692 0 : volume += eta(i, j, k, 0) * dv;
693 0 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
694 0 : Set::Scalar normgrad = grad.lpNorm<2>();
695 0 : Set::Scalar da = normgrad * dv;
696 0 : area += da;
697 :
698 0 : Set::Vector mgrad = Numeric::Gradient(mdot, i, j, k, 0, DX);
699 0 : Set::Scalar mnormgrad = mgrad.lpNorm<2>();
700 0 : Set::Scalar dm = mnormgrad * dv;
701 0 : massflux += dm;
702 :
703 0 : });
704 : }
705 : else {
706 220 : amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
707 : {
708 15560 : volume += eta(i, j, k, 0) * dv;
709 15560 : Set::Vector grad = Numeric::Gradient(eta, i, j, k, 0, DX);
710 15560 : Set::Scalar normgrad = grad.lpNorm<2>();
711 15560 : Set::Scalar da = normgrad * dv;
712 15560 : area += da;
713 15560 : });
714 : }
715 : // time dependent pressure data from experimenta -> p = 0.0954521220950523 * exp(15.289993148880678 * t)
716 220 : }
717 : } // namespace Integrator
718 :
719 :
|
