Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add unified-domain to solve porous flow #1001

Draft
wants to merge 81 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
81 commits
Select commit Hold shift + click to select a range
a26d50e
Fix AW example
tulioricci Mar 11, 2024
2759bef
Remove emissivity from PorousVars
tulioricci Mar 11, 2024
c08165a
Update materials
tulioricci Mar 11, 2024
52cc1a0
Update mechanism
tulioricci Mar 11, 2024
246032f
Fix docs, comments and emissivity for wall_model
tulioricci Mar 11, 2024
857459f
docs
tulioricci Mar 11, 2024
6f0b880
Add test for wall model
tulioricci Mar 11, 2024
b049a4b
Update wall_model
tulioricci Mar 11, 2024
28a6d27
Update materials initializer
tulioricci Mar 11, 2024
8ec893a
Update inviscid flux argument
tulioricci Mar 11, 2024
441b352
Update material and properties when coupling with fluid
tulioricci Mar 11, 2024
18177f2
Update wall model to handle IdealSingleGas
tulioricci Mar 11, 2024
f53ffb1
Add test to test_wallmodel
tulioricci Mar 11, 2024
7236893
Add material property to state creation in boundary conditions
tulioricci Mar 11, 2024
aebbe35
Update gas model to handle IdealSingleGas in porous fluid
tulioricci Mar 11, 2024
6db936f
Update materials initializer to handle IdealSingleGas in porous fluid
tulioricci Mar 11, 2024
3aa5461
docs for the materials
tulioricci Mar 12, 2024
213ed3d
docs for inviscid.py
tulioricci Mar 12, 2024
9a041ba
improve docs in initializer
tulioricci Mar 13, 2024
3656771
Remove dim
tulioricci Mar 14, 2024
4e37629
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 14, 2024
4a1dc95
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 14, 2024
5fc7bb3
Is it necessary to have a specific dt calculation for porous flows?
tulioricci Mar 23, 2024
d9154e0
Merge branch 'main' into tulio/api-updates
tulioricci Mar 24, 2024
f6a529a
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 24, 2024
a766f8d
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 24, 2024
531e278
flake8
tulioricci Mar 25, 2024
c11bc5f
pydocstyle
tulioricci Mar 25, 2024
dcc81e8
Keep Y2 oxidation model for now. Postpone it to Darcy flow.
tulioricci Mar 26, 2024
0f07d4a
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 26, 2024
f863d27
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 26, 2024
94873e3
Improve docs, mypy in carbon_fiber and wall_model
tulioricci Mar 26, 2024
6d4163c
flake8
tulioricci Mar 26, 2024
2505217
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 26, 2024
2412be8
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 26, 2024
fdd7fa8
Use tau instead of progress for Y2 model + docs, mypy
tulioricci Mar 26, 2024
c17f04b
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 26, 2024
fde1088
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 26, 2024
a26e5d5
Organize oxidation models; docs
tulioricci Mar 26, 2024
73b04df
Add velocity to porous flow initializer; fix broken documentation
tulioricci Mar 26, 2024
e33e91c
Optional velocity to porous material
tulioricci Mar 27, 2024
03f47dc
Pass arguments to TACOT degradation
tulioricci Mar 27, 2024
7c1694e
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 27, 2024
7138826
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 27, 2024
1dda90e
Add test for TACOT decomposition; fix mistakes in the model
tulioricci Mar 27, 2024
53b9e1b
Add some extra missing mypy
tulioricci Mar 27, 2024
9fc481b
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Mar 27, 2024
a809035
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Mar 27, 2024
963cc96
Some fixes/updates to work with porous slab
tulioricci Mar 28, 2024
43b52a3
Fix test
tulioricci Mar 28, 2024
1a083c1
Force wall to be positive; add option to return whole set of WVs
tulioricci Apr 2, 2024
4b4b80b
Actually, dont do anything to the wall density but only modify the so…
tulioricci Apr 2, 2024
b760671
Merge branch 'main' into tulio/api-updates
tulioricci Apr 10, 2024
72f5f0f
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Apr 10, 2024
3ac4a7a
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Apr 10, 2024
c7ade2f
Merge branch 'main' into tulio/api-updates
tulioricci Apr 23, 2024
5976817
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci Apr 23, 2024
6e8621f
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci Apr 23, 2024
108d8fd
temp
tulioricci Apr 23, 2024
8917682
Fix to fluxes and BCs and related stuff; Add oxidation models; Modify…
tulioricci May 1, 2024
d2abc3e
Merge branch 'main' into tulio/api-updates
tulioricci May 1, 2024
208aeeb
Merge branch 'tulio/api-updates' into tulio/update_mat_intfce
tulioricci May 1, 2024
b686634
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci May 1, 2024
dc2ece6
Fix float tau; Make alpha explicit in transport
tulioricci May 13, 2024
3134455
Merge branch 'main' into tulio/update_mat_intfce
tulioricci May 13, 2024
fbf0140
Fix test and flake8
tulioricci May 13, 2024
9bbe9b1
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci May 13, 2024
3cf706d
flake8
tulioricci May 13, 2024
c873f70
Merge branch 'main' into tulio/update_mat_intfce
tulioricci May 14, 2024
8b326cc
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci May 14, 2024
3aeac39
Merge branch 'main' into tulio/update_mat_intfce
tulioricci May 15, 2024
16bfdf7
Fix production conflict
tulioricci May 15, 2024
18cf45c
Merge branch 'tulio/update_mat_intfce' into tulio/darcy
tulioricci May 15, 2024
0868c8b
There is no need to have a PR prior to this one.. Incorporate the spe…
tulioricci May 15, 2024
057b2f7
Remove updates to carbon -- that will be added as a new PR
tulioricci May 15, 2024
1999b14
Remove useless changes
tulioricci May 15, 2024
3ce4188
Merge branch 'main' into tulio/darcy
tulioricci May 15, 2024
51b671f
Merge branch 'main' into tulio/darcy
tulioricci Jul 9, 2024
2d5e30b
Merge branch 'main' into tulio/darcy
tulioricci Oct 2, 2024
5a68422
flake8, pylint, mypy
tulioricci Oct 2, 2024
0618d14
mypy; fix np.max in test_wallmodel
tulioricci Oct 3, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 55 additions & 7 deletions mirgecom/boundary.py
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,7 @@
from mirgecom.utils import project_from_base
from mirgecom.viscous import viscous_facial_flux_central, viscous_flux
from mirgecom.inviscid import inviscid_facial_flux_rusanov, inviscid_flux
from mirgecom.gas_model import PorousFlowFluidState


def _ldg_bnd_flux_for_grad(internal_quantity, external_quantity):
Expand Down Expand Up @@ -1090,6 +1091,8 @@ def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
# set the normal momentum to 0
mom_bc = self._slip.momentum_bc(state_minus.momentum_density, nhat)

# this doesn't recompute the internal energy but only substracts the
# kinetic energy from the total energy
energy_bc = (
gas_model.eos.internal_energy(state_minus.cv)
+ 0.5*np.dot(mom_bc, mom_bc)/state_minus.mass_density)
Expand Down Expand Up @@ -1203,11 +1206,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
momentum=free_stream_density*free_stream_velocity,
species_mass=free_stream_spec_mass)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=cv_infinity, gas_model=gas_model,
temperature_seed=free_stream_temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta)
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
"""Return BC fluid state."""
Expand Down Expand Up @@ -1338,11 +1347,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
energy=total_energy,
species_mass=state_minus.cv.species_mass)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=cv_outflow, gas_model=gas_model,
temperature_seed=state_minus.temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta)
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
"""Return state."""
Expand Down Expand Up @@ -1381,11 +1396,17 @@ def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
momentum=state_minus.momentum_density,
species_mass=state_minus.species_mass_density)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=cv_plus, gas_model=gas_model,
temperature_seed=state_minus.temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta)
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def temperature_bc(self, dcoll, dd_bdry, state_minus, **kwargs):
"""Get temperature value used in grad(T)."""
Expand Down Expand Up @@ -1484,11 +1505,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
momentum=rho_boundary * velocity_boundary,
species_mass=species_mass_boundary)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
temperature_seed=state_minus.temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta)
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
"""Return BC fluid state."""
Expand Down Expand Up @@ -1601,11 +1628,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
momentum=rho_boundary*velocity_boundary,
species_mass=species_mass_boundary)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
temperature_seed=state_minus.temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta)
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
"""Return BC fluid state."""
Expand Down Expand Up @@ -1966,11 +1999,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
momentum=mass*velocity,
species_mass=mass*species_mass_fracs)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
temperature_seed=state_minus.temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta)
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
"""Return BC fluid state."""
Expand Down Expand Up @@ -2068,8 +2107,17 @@ def state_plus(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
momentum=mass*velocity,
species_mass=mass*self._spec_mass_fracs)

if isinstance(state_minus, PorousFlowFluidState):
material_densities = state_minus.wv.material_densities
else:
material_densities = None

return make_fluid_state(cv=boundary_cv, gas_model=gas_model,
temperature_seed=state_minus.temperature)
temperature_seed=state_minus.temperature,
smoothness_mu=state_minus.smoothness_mu,
smoothness_kappa=state_minus.smoothness_kappa,
smoothness_beta=state_minus.smoothness_beta,
material_densities=material_densities)

def state_bc(self, dcoll, dd_bdry, gas_model, state_minus, **kwargs):
"""Return BC fluid state."""
Expand Down
57 changes: 38 additions & 19 deletions mirgecom/gas_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -426,7 +426,6 @@ def make_fluid_state(cv, gas_model,
# FIXME per previous review, think of a way to de-couple wall and fluid.
# ~~~ we need to squeeze wall_eos in gas_model because this is easily
# accessible everywhere in the code

tau = gas_model.decomposition_progress(material_densities)
wv = PorousWallVars(
material_densities=material_densities,
Expand All @@ -437,25 +436,45 @@ def make_fluid_state(cv, gas_model,
tortuosity=gas_model.wall_eos.tortuosity(tau=tau)
)

temperature = gas_model.get_temperature(cv=cv, wv=wv,
tseed=temperature_seed)

pressure = gas_model.get_pressure(cv, wv, temperature)
pressure = None
temperature = None

if limiter_func:
cv = limiter_func(cv=cv, wv=wv, pressure=pressure,
temperature=temperature, dd=limiter_dd)
rv = limiter_func(cv=cv, wv=wv, temperature_seed=temperature_seed,
gas_model=gas_model, dd=limiter_dd)
if isinstance(rv, np.ndarray):
cv, pressure, temperature = rv
else:
cv = rv

dv = MixtureDependentVars(
temperature=temperature,
pressure=pressure,
speed_of_sound=gas_model.eos.sound_speed(cv, temperature),
smoothness_mu=smoothness_mu,
smoothness_kappa=smoothness_kappa,
smoothness_d=smoothness_d,
smoothness_beta=smoothness_beta,
species_enthalpies=gas_model.eos.species_enthalpies(cv, temperature),
)
if temperature is None:
temperature = gas_model.get_temperature(cv=cv, wv=wv,
tseed=temperature_seed)
if pressure is None:
pressure = gas_model.get_pressure(cv, wv, temperature)

from mirgecom.eos import MixtureEOS
if isinstance(gas_model.eos, MixtureEOS):
dv = MixtureDependentVars(
temperature=temperature,
pressure=pressure,
speed_of_sound=gas_model.eos.sound_speed(cv, temperature),
smoothness_mu=smoothness_mu,
smoothness_kappa=smoothness_kappa,
smoothness_d=smoothness_d,
smoothness_beta=smoothness_beta,
species_enthalpies=gas_model.eos.species_enthalpies(cv, temperature)
)
else:
dv = GasDependentVars(
temperature=temperature,
pressure=pressure,
speed_of_sound=gas_model.eos.sound_speed(cv, temperature),
smoothness_mu=smoothness_mu,
smoothness_kappa=smoothness_kappa,
smoothness_d=smoothness_d,
smoothness_beta=smoothness_beta
)

tv = gas_model.transport.transport_vars(cv=cv, dv=dv, wv=wv,
eos=gas_model.eos, wall_eos=gas_model.wall_eos)
Expand Down Expand Up @@ -513,7 +532,7 @@ def project_fluid_state(dcoll, src, tgt, state, gas_model, limiter_func=None,
cv_sd = op.project(dcoll, src, tgt, state.cv)

temperature_seed = None
if state.is_mixture:
if state.is_mixture or isinstance(gas_model, PorousFlowModel):
temperature_seed = op.project(dcoll, src, tgt, state.dv.temperature)

if entropy_stable:
Expand Down Expand Up @@ -769,7 +788,7 @@ def make_operator_fluid_states(
]

tseed_interior_pairs = None
if volume_state.is_mixture:
if volume_state.is_mixture or isinstance(gas_model, PorousFlowModel):
# If this is a mixture, we need to exchange the temperature field because
# mixture pressure (used in the inviscid flux calculations) depends on
# temperature and we need to seed the temperature calculation for the
Expand Down
74 changes: 65 additions & 9 deletions mirgecom/inviscid.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,17 +57,27 @@
ConservedVars
)
from mirgecom.utils import normalize_boundaries
from mirgecom.gas_model import PorousFlowFluidState

from arraycontext import outer
from meshmode.dof_array import DOFArray


def inviscid_flux(state):
def inviscid_flux(state, gas_model=None):
r"""Compute the inviscid flux vectors from fluid conserved vars *cv*.

The inviscid fluxes are
$(\rho\vec{V},(\rho{E}+p)\vec{V},\rho(\vec{V}\otimes\vec{V})
+p\mathbf{I}, \rho{Y_s}\vec{V})$
The inviscid fluxes for pure-fluids and porous flows are, respectively,
$$
\begin{bmatrix}(\rho{v})_{j}\\
\left(\rho{E}+p\right){v}_{j}\\
\left((\rho{v})_{j}{v}_{i}+p\delta_{ij}\right)\\
(\rho{Y})_{\alpha}{v}_{j}\end{bmatrix}
\text{ and }
\begin{bmatrix}(\epsilon \rho{v})_{j}\\
\left(\epsilon \rho{E}+\epsilon p\right){v}_{j}\\
\left((\epsilon \rho{v})_{j}{v}_{i}+p\delta_{ij}\right)\\
(\epsilon \rho{Y})_{\alpha}{v}_{j}\end{bmatrix} \text{ .}
$$

.. note::

Expand All @@ -82,6 +92,12 @@ def inviscid_flux(state):

Full fluid conserved and thermal state.

gas_model: :class:`~mirgecom.gas_model.GasModel`

Physical gas model including equation of state, transport,
and kinetic properties as required by fluid state.
Only necessary for porous media flows using the unified-domain.

Returns
-------
:class:`~mirgecom.fluid.ConservedVars`
Expand All @@ -90,11 +106,25 @@ def inviscid_flux(state):
conservation equation.
"""
mass_flux = state.momentum_density
energy_flux = state.velocity * (state.energy_density + state.pressure)

# only the gas part goes into the energy flux
if isinstance(state, PorousFlowFluidState):
gas_energy = state.mass_density*(
gas_model.eos.get_internal_energy(
temperature=state.temperature,
species_mass_fractions=state.cv.species_mass_fractions)
+ 0.5*np.dot(state.velocity, state.velocity)
)
energy_flux = state.velocity * (gas_energy
+ state.wv.void_fraction*state.pressure)
else:
energy_flux = state.velocity * (state.energy_density + state.pressure)

mom_flux = (
state.mass_density * np.outer(state.velocity, state.velocity)
+ np.eye(state.dim)*state.pressure
)

species_mass_flux = \
state.velocity*state.species_mass_density.reshape(-1, 1)

Expand Down Expand Up @@ -146,10 +176,36 @@ def inviscid_facial_flux_rusanov(state_pair, gas_model, normal):
actx = state_pair.int.array_context
lam = actx.np.maximum(state_pair.int.wavespeed, state_pair.ext.wavespeed)
from mirgecom.flux import num_flux_lfr
return num_flux_lfr(f_minus_normal=inviscid_flux(state_pair.int)@normal,
f_plus_normal=inviscid_flux(state_pair.ext)@normal,
q_minus=state_pair.int.cv,
q_plus=state_pair.ext.cv, lam=lam)

# Use only the gas energy in penalization terms
# TODO think of a better way to generalize this for different fluxes??
if isinstance(state_pair.int, PorousFlowFluidState):
gas_energy_int = state_pair.int.mass_density*(
gas_model.eos.get_internal_energy(
temperature=state_pair.int.temperature,
species_mass_fractions=state_pair.int.cv.species_mass_fractions)
+ 0.5*np.dot(state_pair.int.velocity, state_pair.int.velocity)
)
q_minus = state_pair.int.cv.replace(energy=gas_energy_int)

gas_energy_ext = state_pair.ext.mass_density*(
gas_model.eos.get_internal_energy(
temperature=state_pair.ext.temperature,
species_mass_fractions=state_pair.ext.cv.species_mass_fractions)
+ 0.5*np.dot(state_pair.ext.velocity, state_pair.ext.velocity)
)
q_plus = state_pair.ext.cv.replace(energy=gas_energy_ext)

return num_flux_lfr(
f_minus_normal=inviscid_flux(state_pair.int, gas_model=gas_model)@normal,
f_plus_normal=inviscid_flux(state_pair.ext, gas_model=gas_model)@normal,
q_minus=q_minus, q_plus=q_plus, lam=lam)

return num_flux_lfr(
f_minus_normal=inviscid_flux(state_pair.int, gas_model=gas_model)@normal,
f_plus_normal=inviscid_flux(state_pair.ext, gas_model=gas_model)@normal,
q_minus=state_pair.int.cv,
q_plus=state_pair.ext.cv, lam=lam)


def inviscid_facial_flux_hll(state_pair, gas_model, normal):
Expand Down
Loading