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.

Sign up for GitHub

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

Fix diffusion units #177

Closed
wants to merge 3 commits into from
Closed

Fix diffusion units #177

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fb27055
Add check for diffusion units into initialization checks
emmetfrancis Aug 17, 2024
3254d41
Merge branch 'emmetfrancis/param-from-xdmf' into emmetfrancis/fix-dif…
emmetfrancis Aug 17, 2024
2c8bb87
add back D_constant (oops)
emmetfrancis Aug 17, 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
40 changes: 21 additions & 19 deletions smart/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@
)

from .solvers import smartSNESProblem
from .units import unit
from . import unit


logger = logging.getLogger(__name__)
Expand Down Expand Up @@ -632,14 +632,21 @@ def _init_2_4_check_for_unused_parameters_species_compartments(self):
raise ValueError(print_str)

def _init_2_5_link_compartments_to_species(self):
"""Linking compartments and compartment dimensionality to species"""
"""Linking compartments and compartment dimensionality to species,
check for consistency of diffusion coefficient definition"""
logger.debug(
"Linking compartments and compartment dimensionality to species",
extra=dict(format_type="log"),
)
for species in self.sc:
species.compartment = self.cc[species.compartment_name]
species.dimensionality = self.cc[species.compartment_name].dimensionality
# convert diffusion coeff to consistent with mesh
diffusion_conversion = species.diffusion_.to(
species.compartment.compartment_**2 / unit.s
)
species.diffusion_ = species.compartment.compartment_**2 / unit.s
species.D *= diffusion_conversion.magnitude

def _init_2_6_link_species_to_compartments(self):
"""Links species to compartments - a species is considered to be
Expand Down Expand Up @@ -1077,7 +1084,7 @@ def _init_5_2_create_variational_forms(self):
for flux in self.fc:
# -1 factor in flux.form means this is a lhs term
form_type = "boundary_reaction" if flux.is_boundary_condition else "domain_reaction"
flux_form_units = flux.equation_units * flux.measure_units
flux_form_ = flux.equation_ * flux.measure_
# Determine if flux is linear w.r.t. compartment functions
# Use flux.is_linear_wrt_comp and combine with linear_wrt_comp
# (prioritizing former). If compartment is not relevant to flux then it is linear
Expand All @@ -1090,7 +1097,7 @@ def _init_5_2_create_variational_forms(self):
flux.form,
flux.destination_species,
form_type,
flux_form_units,
flux_form_,
True,
linearity_dict,
)
Expand All @@ -1105,15 +1112,15 @@ def _init_5_2_create_variational_forms(self):
v = species.v
D = species.D
dx = species.compartment.mesh.dx
Dform_units = (
species.diffusion_units
* species.concentration_units
* species.compartment.compartment_units ** (species.compartment.dimensionality - 2)
Dform_ = (
species.diffusion_
* species.concentration_
* species.compartment.compartment_ ** (species.compartment.dimensionality - 2)
)
mass_form_units = (
species.concentration_units
mass_form_ = (
species.concentration_
/ unit.s
* species.compartment.compartment_units**species.compartment.dimensionality
* species.compartment.compartment_**species.compartment.dimensionality
)
# diffusion term
if species.D == 0:
Expand All @@ -1123,11 +1130,6 @@ def _init_5_2_create_variational_forms(self):
extra=dict(format_type="log"),
)
else:
if Dform_units != mass_form_units: # unit conversion for consistency
diffusion_conversion = species.diffusion_units.to(
species.compartment.compartment_units**2 / unit.s
)
D *= diffusion_conversion.magnitude
D_constant = d.Constant(D, name=f"D_{species.name}")
if self.config.flags["axisymmetric_model"]:
Dform = x[0] * D_constant * d.inner(d.grad(u), d.grad(v)) * dx
Expand All @@ -1141,7 +1143,7 @@ def _init_5_2_create_variational_forms(self):
Dform,
species,
"diffusion",
Dform_units,
Dform_,
True,
linear_wrt_comp,
)
Expand All @@ -1158,7 +1160,7 @@ def _init_5_2_create_variational_forms(self):
Muform,
species,
"mass_u",
mass_form_units,
mass_form_,
True,
linear_wrt_comp,
)
Expand All @@ -1173,7 +1175,7 @@ def _init_5_2_create_variational_forms(self):
Munform,
species,
"mass_un",
mass_form_units,
mass_form_,
True,
linear_wrt_comp,
)
Expand Down