From cb66f79b4d30356a36348d4fd6ddea8b42e22e9e Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Thu, 12 Jul 2018 16:46:04 +0100 Subject: [PATCH 1/6] Started on reactive problem --- compressible/eos.py | 21 ++++++++++-- compressible/simulation.py | 55 +++++++++++++++++--------------- compressible_react/_defaults | 3 ++ compressible_react/burning.py | 3 ++ compressible_react/simulation.py | 33 +++++++++++++++++-- 5 files changed, 84 insertions(+), 31 deletions(-) create mode 100644 compressible_react/burning.py diff --git a/compressible/eos.py b/compressible/eos.py index 98132e4de..530fb1625 100644 --- a/compressible/eos.py +++ b/compressible/eos.py @@ -24,7 +24,7 @@ def pres(gamma, dens, eint): The pressure """ - p = dens*eint*(gamma - 1.0) + p = dens * eint * (gamma - 1.0) return p @@ -48,7 +48,7 @@ def dens(gamma, pres, eint): The density """ - dens = pres/(eint*(gamma - 1.0)) + dens = pres / (eint * (gamma - 1.0)) return dens @@ -69,5 +69,20 @@ def rhoe(gamma, pres): The internal energy density, rho e """ - rhoe = pres/(gamma - 1.0) + rhoe = pres / (gamma - 1.0) return rhoe + + +def temp(eint, Cv): + """ + Given the internal energy and the heat capacity at constant volume, return the temperature + + Parameters + ---------- + eint : float + internal energy + Cv : float + Heat capacity at constant volume + """ + + return eint / Cv diff --git a/compressible/simulation.py b/compressible/simulation.py index 42169d293..8e80f314d 100644 --- a/compressible/simulation.py +++ b/compressible/simulation.py @@ -20,6 +20,7 @@ class Variables(object): a container class for easy access to the different compressible variable by an integer key """ + def __init__(self, myd): self.nvar = len(myd.names) @@ -58,19 +59,19 @@ def cons_to_prim(U, gamma, ivars, myg): q = myg.scratch_array(nvar=ivars.nq) q[:, :, ivars.irho] = U[:, :, ivars.idens] - q[:, :, ivars.iu] = U[:, :, ivars.ixmom]/U[:, :, ivars.idens] - q[:, :, ivars.iv] = U[:, :, ivars.iymom]/U[:, :, ivars.idens] + q[:, :, ivars.iu] = U[:, :, ivars.ixmom] / U[:, :, ivars.idens] + q[:, :, ivars.iv] = U[:, :, ivars.iymom] / U[:, :, ivars.idens] e = (U[:, :, ivars.iener] - - 0.5*q[:, :, ivars.irho]*(q[:, :, ivars.iu]**2 + - q[:, :, ivars.iv]**2))/q[:, :, ivars.irho] + 0.5 * q[:, :, ivars.irho] * (q[:, :, ivars.iu]**2 + + q[:, :, ivars.iv]**2)) / q[:, :, ivars.irho] q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e) if ivars.naux > 0: - for nq, nu in zip(range(ivars.ix, ivars.ix+ivars.naux), - range(ivars.irhox, ivars.irhox+ivars.naux)): - q[:, :, nq] = U[:, :, nu]/q[:, :, ivars.irho] + for nq, nu in zip(range(ivars.ix, ivars.ix + ivars.naux), + range(ivars.irhox, ivars.irhox + ivars.naux)): + q[:, :, nq] = U[:, :, nu] / q[:, :, ivars.irho] return q @@ -81,18 +82,19 @@ def prim_to_cons(q, gamma, ivars, myg): U = myg.scratch_array(nvar=ivars.nvar) U[:, :, ivars.idens] = q[:, :, ivars.irho] - U[:, :, ivars.ixmom] = q[:, :, ivars.iu]*U[:, :, ivars.idens] - U[:, :, ivars.iymom] = q[:, :, ivars.iv]*U[:, :, ivars.idens] + U[:, :, ivars.ixmom] = q[:, :, ivars.iu] * U[:, :, ivars.idens] + U[:, :, ivars.iymom] = q[:, :, ivars.iv] * U[:, :, ivars.idens] rhoe = eos.rhoe(gamma, q[:, :, ivars.ip]) - U[:, :, ivars.iener] = rhoe + 0.5*q[:, :, ivars.irho]*(q[:, :, ivars.iu]**2 + - q[:, :, ivars.iv]**2) + U[:, :, ivars.iener] = rhoe + \ + 0.5 * q[:, :, ivars.irho] * (q[:, :, ivars.iu]**2 + + q[:, :, ivars.iv]**2) if ivars.naux > 0: - for nq, nu in zip(range(ivars.ix, ivars.ix+ivars.naux), - range(ivars.irhox, ivars.irhox+ivars.naux)): - U[:, :, nu] = q[:, :, nq]*q[:, :, ivars.irho] + for nq, nu in zip(range(ivars.ix, ivars.ix + ivars.naux), + range(ivars.irhox, ivars.irhox + ivars.naux)): + U[:, :, nu] = q[:, :, nq] * q[:, :, ivars.irho] return U @@ -113,7 +115,8 @@ def initialize(self, extra_vars=None, ng=4): # define solver specific boundary condition routines bnd.define_bc("hse", BC.user, is_solid=False) - bnd.define_bc("ramp", BC.user, is_solid=False) # for double mach reflection problem + # for double mach reflection problem + bnd.define_bc("ramp", BC.user, is_solid=False) bc, bc_xodd, bc_yodd = bc_setup(self.rp) @@ -182,10 +185,10 @@ def method_compute_timestep(self): u, v, cs = self.cc_data.get_var(["velocity", "soundspeed"]) # the timestep is min(dx/(|u| + cs), dy/(|v| + cs)) - xtmp = self.cc_data.grid.dx/(abs(u) + cs) - ytmp = self.cc_data.grid.dy/(abs(v) + cs) + xtmp = self.cc_data.grid.dx / (abs(u) + cs) + ytmp = self.cc_data.grid.dy / (abs(v) + cs) - self.dt = cfl*float(min(xtmp.min(), ytmp.min())) + self.dt = cfl * float(min(xtmp.min(), ytmp.min())) def evolve(self): """ @@ -211,19 +214,19 @@ def evolve(self): old_ymom = ymom.copy() # conservative update - dtdx = self.dt/myg.dx - dtdy = self.dt/myg.dy + dtdx = self.dt / myg.dx + dtdy = self.dt / myg.dy for n in range(self.ivars.nvar): var = self.cc_data.get_var_by_index(n) var.v()[:, :] += \ - dtdx*(Flux_x.v(n=n) - Flux_x.ip(1, n=n)) + \ - dtdy*(Flux_y.v(n=n) - Flux_y.jp(1, n=n)) + dtdx * (Flux_x.v(n=n) - Flux_x.ip(1, n=n)) + \ + dtdy * (Flux_y.v(n=n) - Flux_y.jp(1, n=n)) # gravitational source terms - ymom[:, :] += 0.5*self.dt*(dens[:, :] + old_dens[:, :])*grav - ener[:, :] += 0.5*self.dt*(ymom[:, :] + old_ymom[:, :])*grav + ymom[:, :] += 0.5 * self.dt * (dens[:, :] + old_dens[:, :]) * grav + ener[:, :] += 0.5 * self.dt * (ymom[:, :] + old_ymom[:, :]) * grav if self.particles is not None: self.particles.update_particles(self.dt) @@ -257,7 +260,7 @@ def dovis(self): u = q[:, :, ivars.iu] v = q[:, :, ivars.iv] p = q[:, :, ivars.ip] - e = eos.rhoe(gamma, p)/rho + e = eos.rhoe(gamma, p) / rho magvel = np.sqrt(u**2 + v**2) @@ -297,7 +300,7 @@ def dovis(self): # plot particles ax.scatter(particle_positions[:, 0], - particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys") + particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys") ax.set_xlim([myg.xmin, myg.xmax]) ax.set_ylim([myg.ymin, myg.ymax]) diff --git a/compressible_react/_defaults b/compressible_react/_defaults index 543f6d87c..38ef9f702 100644 --- a/compressible_react/_defaults +++ b/compressible_react/_defaults @@ -4,6 +4,7 @@ cfl = 0.8 [eos] gamma = 1.4 ; pres = rho ener (gamma - 1) +cv = 1 ; Heat capacity at constant volume [compressible] @@ -24,3 +25,5 @@ grav = 0.0 ; gravitational acceleration (in y-direction) riemann = HLLC ; HLLC or CGF +[diffusion] +k = 1.0 ; conductivity diff --git a/compressible_react/burning.py b/compressible_react/burning.py new file mode 100644 index 000000000..570e7c2f6 --- /dev/null +++ b/compressible_react/burning.py @@ -0,0 +1,3 @@ +def compute_energy_gen_rate(cc_data, T): + myg = cc_data.grid + return myg.scratch_array() diff --git a/compressible_react/simulation.py b/compressible_react/simulation.py index 8b68c791a..dfabbf314 100644 --- a/compressible_react/simulation.py +++ b/compressible_react/simulation.py @@ -6,6 +6,8 @@ import compressible import compressible.eos as eos +import compressible_react.burning as burning + import util.plot_tools as plot_tools @@ -19,24 +21,51 @@ def initialize(self): """ super().initialize(extra_vars=["fuel", "ash"]) + myd = self.cc_data + def burn(self, dt): """ react fuel to ash """ + + e = self.cc_data.get_var("eint") + ener = self.cc_data.get_var("eint") + dens = self.cc_data.get_var("density") + fuel = self.cc_data.get_var("fuel") + ash = self.cc_data.get_var("ash") + # compute T + Cv = self.rp.get_param("eos.cv") + temp = eos.temp(e, Cv) # compute energy generation rate + omega_dot = burning.compute_energy_gen_rate(self.cc_data, temp) # update energy due to reaction - pass + ener[:,:] += dens * omega_dot * dt + + # update fuel and ash?? def diffuse(self, dt): """ diffuse for dt """ + myg = self.cc_data.grid + + e = self.cc_data.get_var("eint") + ener = self.cc_data.get_var("eint") + # compute T + Cv = self.rp.get_param("eos.cv") + temp = eos.temp(e, Cv) # compute div kappa grad T + k = self.rp.get_param("diffusion.k") + + div_kappa_grad_T = myg.scratch_array() + div_kappa_grad_T.v()[:, :] = k * ( + (temp.ip(1) + temp.ip(-1) - 2.0*temp.v())/myg.dx**2 + + (temp.jp(1) + temp.jp(-1) - 2.0*temp.v())/myg.dy**2) # update energy due to diffusion - pass + ener[:,:] += div_kappa_grad_T * dt def evolve(self): """ From c26c295760dada46cb26a97870f7cdfe8899780d Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Thu, 12 Jul 2018 17:36:30 +0100 Subject: [PATCH 2/6] Ash/ fuel seem to be advecting, added basic burning routines, diffusion looks reasonable --- compressible_react/_defaults | 1 + compressible_react/burning.py | 39 ++++++++++++++++++++++++++++-- compressible_react/simulation.py | 41 ++++++++++++++++++-------------- 3 files changed, 61 insertions(+), 20 deletions(-) diff --git a/compressible_react/_defaults b/compressible_react/_defaults index 38ef9f702..d381e95c4 100644 --- a/compressible_react/_defaults +++ b/compressible_react/_defaults @@ -27,3 +27,4 @@ riemann = HLLC ; HLLC or CGF [diffusion] k = 1.0 ; conductivity +constant_kappa = 1 ; Is the conductivity constant (1) or constant opacity (0) diff --git a/compressible_react/burning.py b/compressible_react/burning.py index 570e7c2f6..d76de94b9 100644 --- a/compressible_react/burning.py +++ b/compressible_react/burning.py @@ -1,3 +1,38 @@ -def compute_energy_gen_rate(cc_data, T): +def energy_and_species_creation(cc_data, T): + """ + Compute the heat release and species create rate + """ myg = cc_data.grid - return myg.scratch_array() + return myg.scratch_array(), myg.scratch_array() + + +def kappa(cc_data, temp, const, constant=1): + """ + Conductivity + + If constant, just returns the constant defined in the params file. + + Otherwise, it uses the formula for conductivity with constant opacity in the Castro/Microphysics library. + + Parameters + ---------- + cc_data : CellCenterData2d + the cell centered data + temp : ArrayIndexer + temperature + const : float + the diffusion constant or opacity + constant : int + Is the conductivity constant (1) or the opacity constant (0) + """ + myg = cc_data.grid + k = myg.scratch_array() + + if constant == 1: + k[:, :] = const + else: + sigma_SB = 5.6704e-5 # Stefan-Boltzmann in cgs + dens = cc_data.get_var("density") + k[:, :] = (16 * sigma_SB * temp**3) / (3 * const * dens) + + return k diff --git a/compressible_react/simulation.py b/compressible_react/simulation.py index dfabbf314..63dea3301 100644 --- a/compressible_react/simulation.py +++ b/compressible_react/simulation.py @@ -21,8 +21,6 @@ def initialize(self): """ super().initialize(extra_vars=["fuel", "ash"]) - myd = self.cc_data - def burn(self, dt): """ react fuel to ash """ @@ -37,12 +35,13 @@ def burn(self, dt): temp = eos.temp(e, Cv) # compute energy generation rate - omega_dot = burning.compute_energy_gen_rate(self.cc_data, temp) + H, omega_dot = burning.energy_and_species_creation(self.cc_data, temp) # update energy due to reaction - ener[:,:] += dens * omega_dot * dt + ener[:, :] += dens * H * dt - # update fuel and ash?? + fuel[:, :] -= omega_dot * dt + ash[:, :] += omega_dot * dt def diffuse(self, dt): """ diffuse for dt """ @@ -57,15 +56,20 @@ def diffuse(self, dt): temp = eos.temp(e, Cv) # compute div kappa grad T - k = self.rp.get_param("diffusion.k") + kappa = self.rp.get_param("diffusion.k") + const_opacity = self.rp.get_param("diffusion.constant_kappa") + k = burning.kappa(self.cc_data, temp, kappa, const_opacity) div_kappa_grad_T = myg.scratch_array() - div_kappa_grad_T.v()[:, :] = k * ( - (temp.ip(1) + temp.ip(-1) - 2.0*temp.v())/myg.dx**2 + - (temp.jp(1) + temp.jp(-1) - 2.0*temp.v())/myg.dy**2) + + div_kappa_grad_T.v()[:, :] = 0.25 * ( + (k.ip(1) * (temp.ip(2) - temp.v()) - + k.ip(-1) * (temp.v() - temp.ip(-2))) / myg.dx**2 + + (k.jp(1) * (temp.jp(2) - temp.v()) - + k.jp(-1) * (temp.v() - temp.jp(-2))) / myg.dy**2) # update energy due to diffusion - ener[:,:] += div_kappa_grad_T * dt + ener[:, :] += div_kappa_grad_T * dt def evolve(self): """ @@ -74,22 +78,22 @@ def evolve(self): """ # we want to do Strang-splitting here - self.burn(self.dt/2) + self.burn(self.dt / 2) - self.diffuse(self.dt/2) + self.diffuse(self.dt / 2) if self.particles is not None: - self.particles.update_particles(self.dt/2) + self.particles.update_particles(self.dt / 2) # note: this will do the time increment and n increment super().evolve() if self.particles is not None: - self.particles.update_particles(self.dt/2) + self.particles.update_particles(self.dt / 2) - self.diffuse(self.dt/2) + self.diffuse(self.dt / 2) - self.burn(self.dt/2) + self.burn(self.dt / 2) def dovis(self): """ @@ -108,13 +112,14 @@ def dovis(self): # outside of a running simulation. gamma = self.cc_data.get_aux("gamma") - q = compressible.cons_to_prim(self.cc_data.data, gamma, ivars, self.cc_data.grid) + q = compressible.cons_to_prim( + self.cc_data.data, gamma, ivars, self.cc_data.grid) rho = q[:, :, ivars.irho] u = q[:, :, ivars.iu] v = q[:, :, ivars.iv] p = q[:, :, ivars.ip] - e = eos.rhoe(gamma, p)/rho + e = eos.rhoe(gamma, p) / rho X = q[:, :, ivars.ix] From 7be6a95406e3855563ebffab06c5f2b869c31b20 Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Mon, 16 Jul 2018 14:59:38 +0100 Subject: [PATCH 3/6] Added powerlaw burning from Microphysics rep --- compressible_react/_defaults | 15 +++- compressible_react/burning.py | 141 ++++++++++++++++++++++++++++++- compressible_react/simulation.py | 29 +++++-- 3 files changed, 171 insertions(+), 14 deletions(-) diff --git a/compressible_react/_defaults b/compressible_react/_defaults index d381e95c4..6caa4749d 100644 --- a/compressible_react/_defaults +++ b/compressible_react/_defaults @@ -5,7 +5,7 @@ cfl = 0.8 [eos] gamma = 1.4 ; pres = rho ener (gamma - 1) cv = 1 ; Heat capacity at constant volume - +cp = 1 [compressible] use_flattening = 1 ; apply flattening at shocks (1) @@ -26,5 +26,16 @@ riemann = HLLC ; HLLC or CGF [diffusion] -k = 1.0 ; conductivity +k = 1e-2 ; conductivity constant_kappa = 1 ; Is the conductivity constant (1) or constant opacity (0) + + + +[network] +network_type = powerlaw +specific_q_burn = 10 +f_act = 1 +t_burn_ref = 10 +rho_burn_ref = 1.5 +rtilde = 1 +nu = 4 diff --git a/compressible_react/burning.py b/compressible_react/burning.py index d76de94b9..d1ea788d1 100644 --- a/compressible_react/burning.py +++ b/compressible_react/burning.py @@ -1,9 +1,142 @@ -def energy_and_species_creation(cc_data, T): +import numpy as np +import compressible.eos as eos +from scipy.integrate import odeint + + +class Network(object): """ - Compute the heat release and species create rate + Class for holding info about the reaction network. """ - myg = cc_data.grid - return myg.scratch_array(), myg.scratch_array() + + def __init__(self, nspec_evolve=0): + """ Constructor """ + + self.nspec = 0 + self.nspec_evolve = nspec_evolve + self.naux = 0 + + self.spec_names = [] + self.A_ion = [] + self.Z_ion = [] + self.E_binding = [] + + self.nrates = 0 + + self.do_constant_volume_burn = False + self.self_heat = False + self.call_eos_in_rhs = True + self.dT_crit = 1e9 + + def initialize(self, rp): + try: + self.Cv = rp.get_param("eos.cv") + except KeyError: + pass + try: + self.Cp = rp.get_param("eos.cp") + except KeyError: + pass + + def finalize(self): + pass + + def enery_gen_rate(self, dydt): + pass + + def temperature_rhs(self, cc_data, E_nuc): + + myg = cc_data.grid + + rhs = myg.scratch_array() + + if self.self_heat: + if self.do_constant_volume_burn: + rhs[:, :] = E_nuc / self.Cv + + else: # constant pressure + rhs[:, :] = E_nuc / self.Cp + + return rhs + + def energy_and_species_creation(self, cc_data): + + myg = cc_data.grid + return myg.scratch_array(), myg.scratch_array(nvar=self.nspec_evolve) + + +class PowerLaw(Network): + """ + powerlaw network. This is a single-step reaction rate. There are only two species, fuel f and ash a which react through the reaction $f + f \rightarrow a + \gamma$. Baryon conservation requires that $A_f = A_a/2$ and charge conservation requires that $Z_f = Z_a/2$. The reaction rate is a powerlaw in temperature. + """ + + def __init__(self): + """ Constructor """ + super().__init__(nspec_evolve=2) + self.nspec = 3 + self.self_heat = True + + def initialize(self, rp): + super().initialize(rp) + + self.spec_names = ["fuel", "ash", "inert"] + + self.A_ion = np.array([2, 4, 8]) + self.Z_ion = np.array([1, 2, 4]) + self.E_binding = np.array( + [0, rp.get_param("network.specific_q_burn"), 0]) + + self.F_act = rp.get_param("network.f_act") + self.T_burn_ref = rp.get_param("network.t_burn_ref") + self.rho_burn_ref = rp.get_param("network.rho_burn_ref") + self.rtilde = rp.get_param("network.rtilde") + self.nu = rp.get_param("network.nu") + + + def enery_gen_rate(self, dydt): + """ + Computers instantaneous energy generation rate + """ + return np.sum(dydt[:, :, :] * + self.A_ion[np.newaxis, + np.newaxis, :self.nspec_evolve] * + self.E_binding[np.newaxis, + np.newaxis, :self.nspec_evolve], + axis=2) + + def energy_and_species_creation(self, cc_data): + """ + RHS of reactions ODE + """ + + myg = cc_data.grid + + xfueltmp = np.maximum(cc_data.get_var("fuel"), 0) + dens = cc_data.get_var("density") + eint = cc_data.get_var("eint") + + T = myg.scratch_array() + T[:, :] = eos.temp(eint, self.Cv) + + rate = myg.scratch_array() + omega_dot = myg.scratch_array(nvar=self.nspec_evolve) + + mask = (T >= self.F_act * self.T_burn_ref) + + rate[mask] = self.rtilde * dens[mask] / self.rho_burn_ref * \ + xfueltmp[mask]**2 * (T[mask] / self.T_burn_ref)**self.nu + + omega_dot[:, :, 0] = -rate + omega_dot[:, :, 1] = rate + + omega_dot[:, :, :] /= \ + self.A_ion[np.newaxis, + np.newaxis, :self.nspec_evolve] + + E_nuc = self.enery_gen_rate(omega_dot) + + # net_T = self.temperature_rhs(cc_data, E_nuc) + + return E_nuc, omega_dot def kappa(cc_data, temp, const, constant=1): diff --git a/compressible_react/simulation.py b/compressible_react/simulation.py index 63dea3301..040153e6a 100644 --- a/compressible_react/simulation.py +++ b/compressible_react/simulation.py @@ -9,6 +9,7 @@ import compressible_react.burning as burning import util.plot_tools as plot_tools +from util import msg class Simulation(compressible.Simulation): @@ -21,27 +22,39 @@ def initialize(self): """ super().initialize(extra_vars=["fuel", "ash"]) + network = self.rp.get_param("network.network_type") + + if network == "powerlaw": + self.network = burning.PowerLaw() + elif network == "null": + self.network = burning.Network(nspec_evolve=2) + else: + msg.fail("ERROR: Do not recognise network %s" % network) + + self.network.initialize(self.rp) + def burn(self, dt): """ react fuel to ash """ - e = self.cc_data.get_var("eint") - ener = self.cc_data.get_var("eint") + # e = self.cc_data.get_var("eint") + ener = self.cc_data.get_var("energy") dens = self.cc_data.get_var("density") fuel = self.cc_data.get_var("fuel") ash = self.cc_data.get_var("ash") # compute T - Cv = self.rp.get_param("eos.cv") - temp = eos.temp(e, Cv) + # Cv = self.rp.get_param("eos.cv") + # temp = eos.temp(e, Cv) # compute energy generation rate - H, omega_dot = burning.energy_and_species_creation(self.cc_data, temp) + H, omega_dot = self.network.energy_and_species_creation(self.cc_data) # update energy due to reaction ener[:, :] += dens * H * dt - fuel[:, :] -= omega_dot * dt - ash[:, :] += omega_dot * dt + fuel[:, :] += omega_dot[:,:,0] * dt * dens + ash[:, :] += omega_dot[:,:,1] * dt * dens + def diffuse(self, dt): """ diffuse for dt """ @@ -49,7 +62,7 @@ def diffuse(self, dt): myg = self.cc_data.grid e = self.cc_data.get_var("eint") - ener = self.cc_data.get_var("eint") + ener = self.cc_data.get_var("energy") # compute T Cv = self.rp.get_param("eos.cv") From 1a2b4b44b6337253475a1080d297728e98f676ca Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Mon, 16 Jul 2018 16:45:09 +0100 Subject: [PATCH 4/6] Added tests --- compressible_react/burning.py | 73 +++++++++++--- compressible_react/problems/test.py | 37 ++++++++ compressible_react/simulation.py | 5 +- compressible_react/tests/test_burning.py | 115 +++++++++++++++++++++++ 4 files changed, 213 insertions(+), 17 deletions(-) create mode 100644 compressible_react/problems/test.py create mode 100644 compressible_react/tests/test_burning.py diff --git a/compressible_react/burning.py b/compressible_react/burning.py index d1ea788d1..e64ee8ec3 100644 --- a/compressible_react/burning.py +++ b/compressible_react/burning.py @@ -1,15 +1,23 @@ import numpy as np import compressible.eos as eos -from scipy.integrate import odeint class Network(object): """ - Class for holding info about the reaction network. + Implements calculations to do with the reaction network. + + This is a base class which implements the null case (no burning). Other networks should inherit from this one. """ def __init__(self, nspec_evolve=0): - """ Constructor """ + """ + Constructor + + Parameters + ---------- + nspec_evolve : int + Number of species to evolve. + """ self.nspec = 0 self.nspec_evolve = nspec_evolve @@ -28,6 +36,9 @@ def __init__(self, nspec_evolve=0): self.dT_crit = 1e9 def initialize(self, rp): + """ + Initialize the class. Will try to find the heat capacities at constant volume and pressure. + """ try: self.Cv = rp.get_param("eos.cv") except KeyError: @@ -37,13 +48,17 @@ def initialize(self, rp): except KeyError: pass - def finalize(self): - pass - - def enery_gen_rate(self, dydt): - pass - def temperature_rhs(self, cc_data, E_nuc): + """ + Calculates the time derivative of the temperature. + + Parameters + ---------- + cc_data : CellCenterData2d + Cell-centered data + E_nuc : Grid2d + Instantaneous energy generation rate. + """ myg = cc_data.grid @@ -59,6 +74,14 @@ def temperature_rhs(self, cc_data, E_nuc): return rhs def energy_and_species_creation(self, cc_data): + """ + Returns the instantaneous energy generation rate and the species creation rate. + + Parameters + ---------- + cc_data : CellCenterData2d + The cell-centered data + """ myg = cc_data.grid return myg.scratch_array(), myg.scratch_array(nvar=self.nspec_evolve) @@ -66,7 +89,12 @@ def energy_and_species_creation(self, cc_data): class PowerLaw(Network): """ - powerlaw network. This is a single-step reaction rate. There are only two species, fuel f and ash a which react through the reaction $f + f \rightarrow a + \gamma$. Baryon conservation requires that $A_f = A_a/2$ and charge conservation requires that $Z_f = Z_a/2$. The reaction rate is a powerlaw in temperature. + powerlaw network. This is a single-step reaction rate. + There are only two species, fuel f and ash a which react + through the reaction $f + f \rightarrow a + \gamma$. + Baryon conservation requires that $A_f = A_a/2$ and + charge conservation requires that $Z_f = Z_a/2$. The + reaction rate is a powerlaw in temperature. """ def __init__(self): @@ -76,6 +104,14 @@ def __init__(self): self.self_heat = True def initialize(self, rp): + """ + Initialize the object, loading up a number of runtime parameters. + + Parameters + ---------- + rp : RuntimeParameters + Runtime paramters + """ super().initialize(rp) self.spec_names = ["fuel", "ash", "inert"] @@ -91,10 +127,14 @@ def initialize(self, rp): self.rtilde = rp.get_param("network.rtilde") self.nu = rp.get_param("network.nu") - def enery_gen_rate(self, dydt): """ - Computers instantaneous energy generation rate + Computes instantaneous energy generation rate + + Parameters + ---------- + dydt : Grid2d + Species creation rate """ return np.sum(dydt[:, :, :] * self.A_ion[np.newaxis, @@ -105,14 +145,19 @@ def enery_gen_rate(self, dydt): def energy_and_species_creation(self, cc_data): """ - RHS of reactions ODE + Returns the instantaneous energy generation rate and the species creation rate. + + Parameters + ---------- + cc_data : CellCenterData2d + The cell-centered data """ myg = cc_data.grid - xfueltmp = np.maximum(cc_data.get_var("fuel"), 0) dens = cc_data.get_var("density") eint = cc_data.get_var("eint") + xfueltmp = np.maximum(cc_data.get_var("fuel")/dens, 0) T = myg.scratch_array() T[:, :] = eos.temp(eint, self.Cv) diff --git a/compressible_react/problems/test.py b/compressible_react/problems/test.py new file mode 100644 index 000000000..dedd15101 --- /dev/null +++ b/compressible_react/problems/test.py @@ -0,0 +1,37 @@ +from __future__ import print_function + +import sys +import mesh.patch as patch + + +def init_data(my_data, rp): + """ an init routine for unit testing """ + + # make sure that we are passed a valid patch object + if not isinstance(my_data, patch.CellCenterData2d): + print("ERROR: patch invalid in sedov.py") + print(my_data.__class__) + sys.exit() + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + fuel = my_data.get_var("fuel") + ash = my_data.get_var("ash") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + dens[:, :] = 1.0 + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + ener[:, :] = 2.5 + fuel[:, :] = 0.2 + ash[:, :] = 0.8 + + +def finalize(): + """ print out any information to the user at the end of the run """ + pass diff --git a/compressible_react/simulation.py b/compressible_react/simulation.py index 040153e6a..132218d7b 100644 --- a/compressible_react/simulation.py +++ b/compressible_react/simulation.py @@ -52,9 +52,8 @@ def burn(self, dt): # update energy due to reaction ener[:, :] += dens * H * dt - fuel[:, :] += omega_dot[:,:,0] * dt * dens - ash[:, :] += omega_dot[:,:,1] * dt * dens - + fuel[:, :] += omega_dot[:, :, 0] * dt * dens + ash[:, :] += omega_dot[:, :, 1] * dt * dens def diffuse(self, dt): """ diffuse for dt """ diff --git a/compressible_react/tests/test_burning.py b/compressible_react/tests/test_burning.py new file mode 100644 index 000000000..b5644c988 --- /dev/null +++ b/compressible_react/tests/test_burning.py @@ -0,0 +1,115 @@ +import compressible_react.burning as burning +from util import runparams +from numpy.testing import assert_array_equal, assert_array_almost_equal +import numpy as np +import compressible_react.simulation as sn + + +class TestSimulation(object): + @classmethod + def setup_class(cls): + """ this is run once for each class before any tests """ + pass + + @classmethod + def teardown_class(cls): + """ this is run once for each class after all tests """ + pass + + def setup_method(self): + """ this is run before each test """ + self.rp = runparams.RuntimeParameters() + self.rp.params["eos.cv"] = 0.5 + self.rp.params["eos.cp"] = 1 + self.rp.params["eos.gamma"] = 1.4 + self.rp.params["compressible.grav"] = 1.0 + self.rp.params["mesh.nx"] = 8 + self.rp.params["mesh.ny"] = 8 + self.rp.params["particles.do_particles"] = 0 + + self.sim = sn.Simulation("compressible_react", "test", self.rp) + + self.rp = self.sim.rp + + def teardown_method(self): + """ this is run after each test """ + self.rp = None + self.sim = None + + def test_null_network(self): + """ + Test the null network + """ + + self.rp.params["network.network_type"] = "null" + self.sim.initialize() + + ntk = self.sim.network + + assert ntk.nspec == 0 + assert ntk.spec_names == [] + + ntk.initialize(self.rp) + + assert ntk.Cv == self.rp.params["eos.cv"] + assert ntk.Cp == self.rp.params["eos.cp"] + + H, omega_dot = ntk.energy_and_species_creation(self.sim.cc_data) + + assert_array_equal(H, np.zeros_like(H)) + assert_array_equal(omega_dot, np.zeros_like(omega_dot)) + + def test_powerlaw_network(self): + """ + Test the powerlaw network + """ + + self.rp.params["network.network_type"] = "powerlaw" + self.rp.params["network.specific_q_burn"] = 10 + self.rp.params["network.f_act"] = 1 + self.rp.params["network.t_burn_ref"] = 2 + self.rp.params["network.rho_burn_ref"] = 3 + self.rp.params["network.rtilde"] = 4 + self.rp.params["network.nu"] = 5 + + self.sim.initialize() + + ntk = self.sim.network + + assert ntk.nspec == 3 + assert ntk.nspec_evolve == 2 + assert_array_equal(ntk.A_ion, [2, 4, 8]) + + rate = 125 / 24 + + myg = self.sim.cc_data.grid + omega_dot_correct = myg.scratch_array(nvar=2) + omega_dot_correct[:, :, 0] = -rate / 2 + omega_dot_correct[:, :, 1] = rate / 4 + + E_nuc_correct = myg.scratch_array() + E_nuc_correct[:, :] = rate * self.rp.params["network.specific_q_burn"] + + E_nuc, omega_dot = ntk.energy_and_species_creation(self.sim.cc_data) + + assert_array_almost_equal(omega_dot, omega_dot_correct) + assert_array_almost_equal(E_nuc, E_nuc_correct) + + def test_kappa(self): + """ + Test the conductivity + """ + self.rp.params["network.network_type"] = "null" + self.sim.initialize() + + myd = self.sim.cc_data + myg = myd.grid + + k = burning.kappa(myd, 1, 30, 1) + + assert_array_equal(myg.scratch_array() + 30, k) + + temp = myg.scratch_array() + 1 + k_correct = myg.scratch_array() + 1.512106667e-4 + + assert_array_almost_equal(k_correct, burning.kappa(myd, temp, 2, 0)) From 417cf61780a0f198b0e555e70cab8ca71b71655c Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Tue, 17 Jul 2018 15:50:25 +0100 Subject: [PATCH 5/6] Added docs --- compressible/simulation.py | 9 +- compressible_react/burning.py | 38 ++++--- compressible_react/problems/flame.py | 26 +++-- compressible_react/problems/inputs.flame | 5 + compressible_react/simulation.py | 28 ++++-- docs/source/compressible_react.rst | 6 ++ docs/source/compressible_react_basics.rst | 116 ++++++++++++++++++++++ docs/source/design.rst | 13 ++- docs/source/index.rst | 1 + 9 files changed, 206 insertions(+), 36 deletions(-) create mode 100644 docs/source/compressible_react_basics.rst diff --git a/compressible/simulation.py b/compressible/simulation.py index 8e80f314d..373a9e59b 100644 --- a/compressible/simulation.py +++ b/compressible/simulation.py @@ -207,8 +207,10 @@ def evolve(self): myg = self.cc_data.grid - Flux_x, Flux_y = flx.unsplit_fluxes(self.cc_data, self.aux_data, self.rp, - self.ivars, self.solid, self.tc, self.dt) + Flux_x, Flux_y = flx.unsplit_fluxes(self.cc_data, + self.aux_data, self.rp, + self.ivars, self.solid, + self.tc, self.dt) old_dens = dens.copy() old_ymom = ymom.copy() @@ -300,7 +302,8 @@ def dovis(self): # plot particles ax.scatter(particle_positions[:, 0], - particle_positions[:, 1], s=5, c=colors, alpha=0.8, cmap="Greys") + particle_positions[:, 1], s=5, c=colors, + alpha=0.8, cmap="Greys") ax.set_xlim([myg.xmin, myg.xmax]) ax.set_ylim([myg.ymin, myg.ymax]) diff --git a/compressible_react/burning.py b/compressible_react/burning.py index e64ee8ec3..870898715 100644 --- a/compressible_react/burning.py +++ b/compressible_react/burning.py @@ -23,13 +23,11 @@ def __init__(self, nspec_evolve=0): self.nspec_evolve = nspec_evolve self.naux = 0 - self.spec_names = [] + self.spec_names = None self.A_ion = [] self.Z_ion = [] self.E_binding = [] - self.nrates = 0 - self.do_constant_volume_burn = False self.self_heat = False self.call_eos_in_rhs = True @@ -88,12 +86,12 @@ def energy_and_species_creation(self, cc_data): class PowerLaw(Network): - """ + r""" powerlaw network. This is a single-step reaction rate. - There are only two species, fuel f and ash a which react - through the reaction $f + f \rightarrow a + \gamma$. - Baryon conservation requires that $A_f = A_a/2$ and - charge conservation requires that $Z_f = Z_a/2$. The + There are only two species, fuel :math:`f` and ash :math:`a` which react + through the reaction :math:`f + f \rightarrow a + \gamma`. + Baryon conservation requires that :math:`A_f = A_a/2` and + charge conservation requires that :math:`Z_f = Z_a/2`. The reaction rate is a powerlaw in temperature. """ @@ -128,8 +126,14 @@ def initialize(self, rp): self.nu = rp.get_param("network.nu") def enery_gen_rate(self, dydt): - """ - Computes instantaneous energy generation rate + r""" + Computes instantaneous energy generation rate. It is given by + + .. math:: + + f(\dot{Y}_k) = \sum_k \dot{Y}_k A_k E_{\text{bin},k}, + + where :math:`E_{\text{bin},k}` is the binding energy of species :math:`k`. Parameters ---------- @@ -144,9 +148,19 @@ def enery_gen_rate(self, dydt): axis=2) def energy_and_species_creation(self, cc_data): - """ + r""" Returns the instantaneous energy generation rate and the species creation rate. + The rate is given by + + .. math:: + + \dot{\omega}_k = \tilde{r} \frac{\rho}{\rho_{\text{ref}}}X_k^2 \left(\frac{T}{T_{\text{ref}}}\right)^\nu, + + where :math:`\tilde{r}` is the coefficient for the reaction rate, :math:`\rho_{\text{ref}}` and :math:`T_{\text{ref}}` are the reference density and temperature, and :math:`\nu` is the exponent for the temperature. + + :math:`\tilde{r}` is zero if the temperature is below some activation temperature, given by some fraction :math:`F_{\text{act}}` of the reference temperature. + Parameters ---------- cc_data : CellCenterData2d @@ -184,7 +198,7 @@ def energy_and_species_creation(self, cc_data): return E_nuc, omega_dot -def kappa(cc_data, temp, const, constant=1): +def k_th(cc_data, temp, const, constant=1): """ Conductivity diff --git a/compressible_react/problems/flame.py b/compressible_react/problems/flame.py index 75ad00580..886d5c4ea 100644 --- a/compressible_react/problems/flame.py +++ b/compressible_react/problems/flame.py @@ -23,6 +23,8 @@ def init_data(my_data, rp): xmom = my_data.get_var("x-momentum") ymom = my_data.get_var("y-momentum") ener = my_data.get_var("energy") + fuel = my_data.get_var("fuel") + ash = my_data.get_var("ash") # initialize the components, remember, that ener here is rho*eint # + 0.5*rho*v**2, where eint is the specific internal energy @@ -30,6 +32,8 @@ def init_data(my_data, rp): dens[:, :] = 1.0 xmom[:, :] = 0.0 ymom[:, :] = 0.0 + fuel[:, :] = dens + ash[:, :] = 0.0 E_sedov = 1.0 @@ -44,8 +48,8 @@ def init_data(my_data, rp): ymin = rp.get_param("mesh.ymin") ymax = rp.get_param("mesh.ymax") - xctr = 0.5*(xmin + xmax) - yctr = 0.5*(ymin + ymax) + xctr = 0.5 * (xmin + xmax) + yctr = 0.5 * (ymin + ymax) # initialize the pressure by putting the explosion energy into a # volume of constant pressure. Then compute the energy in a zone @@ -56,31 +60,35 @@ def init_data(my_data, rp): (my_data.grid.y2d - yctr)**2) p = 1.e-5 - ener[:, :] = p/(gamma - 1.0) + ener[:, :] = p / (gamma - 1.0) - for i, j in np.transpose(np.nonzero(dist < 2.0*r_init)): + for i, j in np.transpose(np.nonzero(dist < 2.0 * r_init)): pzone = 0.0 for ii in range(nsub): for jj in range(nsub): - xsub = my_data.grid.xl[i] + (my_data.grid.dx/nsub)*(ii + 0.5) - ysub = my_data.grid.yl[j] + (my_data.grid.dy/nsub)*(jj + 0.5) + xsub = my_data.grid.xl[i] + \ + (my_data.grid.dx / nsub) * (ii + 0.5) + ysub = my_data.grid.yl[j] + \ + (my_data.grid.dy / nsub) * (jj + 0.5) dist = np.sqrt((xsub - xctr)**2 + (ysub - yctr)**2) if dist <= r_init: - p = (gamma - 1.0)*E_sedov/(pi*r_init*r_init) + p = (gamma - 1.0) * E_sedov / (pi * r_init * r_init) else: p = 1.e-5 pzone += p - p = pzone/(nsub*nsub) + p = pzone / (nsub * nsub) - ener[i, j] = p/(gamma - 1.0) + ener[i, j] = p / (gamma - 1.0) + ash[i, j] = dens[i, j] + fuel[i, j] = 0.0 def finalize(): diff --git a/compressible_react/problems/inputs.flame b/compressible_react/problems/inputs.flame index 1dabca855..a1b76354d 100644 --- a/compressible_react/problems/inputs.flame +++ b/compressible_react/problems/inputs.flame @@ -36,3 +36,8 @@ r_init = 0.01 [vis] dovis = 1 + + +[network] +rho_burn_ref = 1 +t_burn_ref = 100 diff --git a/compressible_react/simulation.py b/compressible_react/simulation.py index 132218d7b..0303abbeb 100644 --- a/compressible_react/simulation.py +++ b/compressible_react/simulation.py @@ -20,7 +20,6 @@ def initialize(self): the data is the same as the compressible solver, but we supply additional variables. """ - super().initialize(extra_vars=["fuel", "ash"]) network = self.rp.get_param("network.network_type") @@ -33,19 +32,16 @@ def initialize(self): self.network.initialize(self.rp) + super().initialize(extra_vars=self.network.spec_names) + def burn(self, dt): """ react fuel to ash """ - # e = self.cc_data.get_var("eint") ener = self.cc_data.get_var("energy") dens = self.cc_data.get_var("density") fuel = self.cc_data.get_var("fuel") ash = self.cc_data.get_var("ash") - # compute T - # Cv = self.rp.get_param("eos.cv") - # temp = eos.temp(e, Cv) - # compute energy generation rate H, omega_dot = self.network.energy_and_species_creation(self.cc_data) @@ -55,6 +51,18 @@ def burn(self, dt): fuel[:, :] += omega_dot[:, :, 0] * dt * dens ash[:, :] += omega_dot[:, :, 1] * dt * dens + # enforce 0 <= X_k <= 1 + fuel[fuel < 0] = 0 + ash[ash < 0] = 0 + + fuel[fuel > dens] = dens[fuel > dens] + ash[ash > dens] = dens[ash > dens] + + # enforce sum_k X_k = 1 by renormalizing + sum = (fuel + ash) / dens + fuel[:, :] /= sum + ash[:, :] /= sum + def diffuse(self, dt): """ diffuse for dt """ @@ -68,9 +76,9 @@ def diffuse(self, dt): temp = eos.temp(e, Cv) # compute div kappa grad T - kappa = self.rp.get_param("diffusion.k") + k_const = self.rp.get_param("diffusion.k") const_opacity = self.rp.get_param("diffusion.constant_kappa") - k = burning.kappa(self.cc_data, temp, kappa, const_opacity) + k = burning.k_th(self.cc_data, temp, k_const, const_opacity) div_kappa_grad_T = myg.scratch_array() @@ -139,8 +147,8 @@ def dovis(self): myg = self.cc_data.grid - fields = [rho, magvel, p, e, X] - field_names = [r"$\rho$", r"U", "p", "e", r"$X_\mathrm{fuel}$"] + fields = [rho, magvel, u, p, e, X] + field_names = [r"$\rho$", r"U", "u", "p", "e", r"$X_\mathrm{fuel}$"] f, axes, cbar_title = plot_tools.setup_axes(myg, len(fields)) diff --git a/docs/source/compressible_react.rst b/docs/source/compressible_react.rst index 826bde4f6..bd2d11941 100644 --- a/docs/source/compressible_react.rst +++ b/docs/source/compressible_react.rst @@ -24,4 +24,10 @@ compressible\_react\.simulation module :undoc-members: :show-inheritance: +compressible\_react\.burning module +-------------------------------------- +.. automodule:: compressible_react.burning + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/compressible_react_basics.rst b/docs/source/compressible_react_basics.rst new file mode 100644 index 000000000..74439b5f7 --- /dev/null +++ b/docs/source/compressible_react_basics.rst @@ -0,0 +1,116 @@ +Reactive flow +============= + +The ``compressible_react`` solver can be used for modelling compressible hydrodynamics with reactions. + +The equations of compressible hydrodynamics with reactive and diffusive source terms take the form: + +.. math:: + + \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho U) &= 0 \\ + \frac{\partial \left(\rho X_k\right)}{\partial t} + \nabla \cdot (\rho U X_k) &= \rho \dot{\omega}_k \\ + \frac{\partial (\rho U)}{\partial t} + \nabla \cdot (\rho U U) + \nabla p &= \rho g \\ + \frac{\partial (\rho E)}{\partial t} + \nabla \cdot [(\rho E + p ) U] &= \nabla \cdot \left(k_{\text{th}} \nabla T\right) + \rho H_{\text{nuc}} + +where :math:`X_k = \rho_k/\rho` is the mass fraction of species :math:`k` (and :math:`\sum_k X_k = 1`), :math:`\dot{\omega}_k` is the species creation rate, :math:`k_{\text{th}}` is the thermal conductivity, :math:`T` is the temperature and :math:`H_{\text{nuc}}` is the time-rate of energy release per unit mass. + +Burning and diffusion are implemented using *Strang splitting*. We evolve the burning and diffusion terms by half a timestep, evolve the compressible hydrodynamics terms by a full timestep, then evolve the diffusion and burning terms by the remaining half timestep. This approach makes the operator splitting second-order in time. + + +Burning +------- +In order include reactions in our simulation, we need to define a reaction network. This is handled by the :func:`Network ` class. This class allows you to define the number of species in the network, the species' properties (atomic mass, proton number and binding energy), and how they react. Once the network has been initialized at the beginning of the simulation, the instantaneous energy generation rate and the species creation rate can be found by calling the member function :func:`energy_and_species_creation `. + +The base ``Network`` class implements a 'null' network - no burning or species. A single-step reaction network is implemented in the :func:`PowerLaw ` class. This models a network made up of two species: fuel, :math:`f`, and ash, :math:`a`. They react +through the reaction + +.. math:: + f + f \rightarrow a + \gamma. + +Baryon conservation requires that :math:`A_f = A_a/2` and +charge conservation requires that :math:`Z_f = Z_a/2`. The +reaction rate is a powerlaw in temperature: + +.. math:: + + \dot{\omega}_k = \tilde{r} \frac{\rho}{\rho_{\text{ref}}}X_k^2 \left(\frac{T}{T_{\text{ref}}}\right)^\nu, + +where :math:`\tilde{r}` is the coefficient for the reaction rate, :math:`\rho_{\text{ref}}` and :math:`T_{\text{ref}}` are the reference density and temperature, and :math:`\nu` is the exponent for the temperature. + +:math:`\tilde{r}` is zero if the temperature is below some activation temperature, given by some fraction :math:`F_{\text{act}}` of the reference temperature. + +The reaction network can be configured using the following runtime parameters: + ++--------------------------------------------------------------------------------+ +|``[network]`` | ++=======================+========================================================+ +|``network_type`` | What type of network? ``null`` or ``powerlaw`` | ++-----------------------+--------------------------------------------------------+ +|``f_act`` | activation temperature factor | ++-----------------------+--------------------------------------------------------+ +|``t_burn_ref`` | reference temperature | ++-----------------------+--------------------------------------------------------+ +|``rho_burn_ref`` | reference density | ++-----------------------+--------------------------------------------------------+ +|``rtilde`` | coefficient for the reaction rate | ++-----------------------+--------------------------------------------------------+ +|``nu`` | exponent for the temperature | ++-----------------------+--------------------------------------------------------+ + + + +Diffusion +--------- + +Diffusion follows Fick's law: the heat flux is proportional to the gradient of the temperature, :math:`F_{\text{cond}} = -k_{\text{th}}\nabla T`. The thermal conductivity :math:`-k_{\text{th}}` can be calculated in a number of ways. The simplest is to set it to be constant everywhere. Alternatively, we can set the opacity :math:`\kappa` to be constant. In this case, the conductivity is given by + +.. math:: + + k_{\text{th}} = \frac{16 \sigma_{\text{SB}} T^3}{3 \kappa \rho}, + +where :math:`\sigma_{\text{SB}}` is the Stefan-Boltzmann constant. + +The diffusion can be configured using the following runtime parameters: + ++--------------------------------------------------------------------------------+ +|``[diffusion]`` | ++=======================+========================================================+ +|``k`` | conductivity constant | ++-----------------------+--------------------------------------------------------+ +|``constant_kappa`` | Constant conductivity (1) or opacity (0)? | ++-----------------------+--------------------------------------------------------+ + + +Code examples +------------- + +To include reactions in the code, we need to first create and initialize the reaction network. We do this by creating a ``Network`` object (or an object of a subclass), then calling the object's :func:`initialize ` function. This is done in the ``Simulation.initialize`` function in ``compressible_react``. + +.. code-block:: python + + from util import runparams + import compressible_react.burning as burning + + # create and import runtime parameters from file + rp = runparams.RuntimeParameters() + + .... + + network = burning.PowerLaw() + network.initialize(rp) + +The internal energy and conserved species mass fractions are updated by calling the :func:`burn ` function. This in turn calls the network object's :func:`energy_and_species_creation ` function to calculate the instantaneous energy generation rate and species creation rates: + +.. code-block:: python + + H, omega_dot = network.energy_and_species_creation(cc_data) + +where ``cc_data`` is a ``CellCenterData2d`` object. + +The diffusion term in the energy equation is found using :func:`diffusion `. This calls the :func:`k_th ` function to calculate the conductivity: + +.. code-block:: python + + k = burning.k_th(cc_data, temp, k_const), + +then computes :math:`\nabla\cdot (k\nabla T)`. diff --git a/docs/source/design.rst b/docs/source/design.rst index 44ff46c49..555b6b586 100644 --- a/docs/source/design.rst +++ b/docs/source/design.rst @@ -54,11 +54,16 @@ The overall structure is: * ``problems/``: The problem setups for the compressible hydrodynamics solver. * ``tests/``: Reference compressible hydro output for regression testing. -* ``compressible_rk/``: The compressible hydrodynamics solver using method of lines integration. +* ``compressible_react/``: The compressible hydrodynamics solver with reactions. - * ``problems/``: This is a symbolic link to the compressible/problems/ directory. + * ``problems/``: The problem setups for the reactive compressible hydrodynamics solver. * ``tests/``: Reference compressible hydro output for regression testing. +* ``compressible_rk/``: The compressible hydrodynamics solver using method of lines integration. + + * ``problems/``: This is a symbolic link to the compressible/problems/ directory. + * ``tests/``: Reference compressible hydro output for regression testing. + * ``diffusion/``: The implicit (thermal) diffusion solver. All diffusion-specific routines live here. * ``problems/``: The problem setups for the diffusion solver. @@ -85,6 +90,10 @@ The overall structure is: * ``problems/``: The problem setups for when the multigrid solver is used in a stand-alone fashion. * ``tests/``: Reference multigrid solver solutions (from when the multigrid solver is used stand-alone) for regression testing. +* ``particles/``: The particles solver + + * ``tests/``: Unit tests for particles solver + * ``swe/``: The shallow water solver. * ``problems/``: The problem setups for the shallow water solver. diff --git a/docs/source/index.rst b/docs/source/index.rst index d2d2e84e6..d87f3a80f 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -38,6 +38,7 @@ pyro: a python hydro code incompressible_basics lowmach_basics particles_basics + compressible_react_basics swe_basics .. toctree:: From 8659010064e42f7030a967c532e7aecd09a81178 Mon Sep 17 00:00:00 2001 From: Alice Harpole Date: Tue, 17 Jul 2018 16:31:48 +0100 Subject: [PATCH 6/6] Fixed tests --- compressible_react/burning.py | 11 ++++++++--- compressible_react/problems/test.py | 14 +++++++++----- compressible_react/tests/test_burning.py | 8 ++++---- 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/compressible_react/burning.py b/compressible_react/burning.py index 870898715..92ef43e9f 100644 --- a/compressible_react/burning.py +++ b/compressible_react/burning.py @@ -157,9 +157,14 @@ def energy_and_species_creation(self, cc_data): \dot{\omega}_k = \tilde{r} \frac{\rho}{\rho_{\text{ref}}}X_k^2 \left(\frac{T}{T_{\text{ref}}}\right)^\nu, - where :math:`\tilde{r}` is the coefficient for the reaction rate, :math:`\rho_{\text{ref}}` and :math:`T_{\text{ref}}` are the reference density and temperature, and :math:`\nu` is the exponent for the temperature. + where :math:`\tilde{r}` is the coefficient for the reaction rate, + :math:`\rho_{\text{ref}}` and :math:`T_{\text{ref}}` are the reference + density and temperature, and :math:`\nu` is the exponent for the + temperature. - :math:`\tilde{r}` is zero if the temperature is below some activation temperature, given by some fraction :math:`F_{\text{act}}` of the reference temperature. + :math:`\tilde{r}` is zero if the temperature is below some activation + temperature, given by some fraction :math:`F_{\text{act}}` of the + reference temperature. Parameters ---------- @@ -171,7 +176,7 @@ def energy_and_species_creation(self, cc_data): dens = cc_data.get_var("density") eint = cc_data.get_var("eint") - xfueltmp = np.maximum(cc_data.get_var("fuel")/dens, 0) + xfueltmp = np.maximum(cc_data.get_var("fuel") / dens, 0) T = myg.scratch_array() T[:, :] = eos.temp(eint, self.Cv) diff --git a/compressible_react/problems/test.py b/compressible_react/problems/test.py index dedd15101..51e4be9e8 100644 --- a/compressible_react/problems/test.py +++ b/compressible_react/problems/test.py @@ -18,9 +18,6 @@ def init_data(my_data, rp): xmom = my_data.get_var("x-momentum") ymom = my_data.get_var("y-momentum") ener = my_data.get_var("energy") - fuel = my_data.get_var("fuel") - ash = my_data.get_var("ash") - # initialize the components, remember, that ener here is rho*eint # + 0.5*rho*v**2, where eint is the specific internal energy # (erg/g) @@ -28,8 +25,15 @@ def init_data(my_data, rp): xmom[:, :] = 0.0 ymom[:, :] = 0.0 ener[:, :] = 2.5 - fuel[:, :] = 0.2 - ash[:, :] = 0.8 + + try: + fuel = my_data.get_var("fuel") + ash = my_data.get_var("ash") + + fuel[:, :] = 0.2 + ash[:, :] = 0.8 + except IndexError: + pass def finalize(): diff --git a/compressible_react/tests/test_burning.py b/compressible_react/tests/test_burning.py index b5644c988..cedd80cc2 100644 --- a/compressible_react/tests/test_burning.py +++ b/compressible_react/tests/test_burning.py @@ -47,7 +47,7 @@ def test_null_network(self): ntk = self.sim.network assert ntk.nspec == 0 - assert ntk.spec_names == [] + assert ntk.spec_names is None ntk.initialize(self.rp) @@ -95,7 +95,7 @@ def test_powerlaw_network(self): assert_array_almost_equal(omega_dot, omega_dot_correct) assert_array_almost_equal(E_nuc, E_nuc_correct) - def test_kappa(self): + def test_k_th(self): """ Test the conductivity """ @@ -105,11 +105,11 @@ def test_kappa(self): myd = self.sim.cc_data myg = myd.grid - k = burning.kappa(myd, 1, 30, 1) + k = burning.k_th(myd, 1, 30, 1) assert_array_equal(myg.scratch_array() + 30, k) temp = myg.scratch_array() + 1 k_correct = myg.scratch_array() + 1.512106667e-4 - assert_array_almost_equal(k_correct, burning.kappa(myd, temp, 2, 0)) + assert_array_almost_equal(k_correct, burning.k_th(myd, temp, 2, 0))