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