From 8aa4f3e69387cbfe2adc956ca5a548783706a8ac Mon Sep 17 00:00:00 2001 From: Jim Schaff Date: Sun, 27 Sep 2026 01:33:52 -0400 Subject: [PATCH] Region variables: model VolumeRegion/MembraneRegion equations, read, write, and results VCell writes a well-mixed species in a spatial application, and a whole-membrane quantity such as the membrane potential, as a region variable: one value per connected region, with a VolumeRegionEquation or MembraneRegionEquation. pyvcell dropped those equations on read, so consumers saw undeclared variables. Result processing also crashed on any model with one, because a volume function referencing it had no numexpr binding (KeyError). - models_math: VolumeRegionEquation (uniform_rate, volume_rate, initial, solution_type, solution) on CompartmentSubDomain, and MembraneRegionEquation (with membrane_rate) on MembraneSubDomain. Both are exported from pyvcell.vcml. - vcml_reader / vcml_writer: read and write them as VCell's XmlReader and Xmlproducer do: UniformRate and the region rate (default 0.0), Initial, and , whose text is an exact solution. - SpeciesMapping: force_constant, well_mixed and force_continuous, the LocalizedCompoundSpec flags. The writer used to drop them, so a round trip (or vc.simulate) turned a well-mixed species spatial and unclamped a clamped species. - zarr_writer: each domain-qualified VOLUME_REGION state variable is replicated onto every voxel of its region, values[volume_region_map]. It becomes a channel after the volume variables and a binding for the functions, as x/y/z/t already are. Membrane-region variables stay out of the zarr, like membrane variables (it is a volume grid). Tests use a new 20 KB fixture: SmallSpatialProject_3D with s1 (volume) and s2 (membrane) well-mixed, math regenerated by libvcell. They cover the region equations as read, the math and flag round trips, and an fvsolver run whose s1 channel is constant across subdomain0 and equals the solver's per-region value at every output time. Co-Authored-By: Claude Opus 5.5 (1M context) --- pyvcell/_internal/simdata/zarr_writer.py | 34 +- pyvcell/vcml/__init__.py | 4 + pyvcell/vcml/models_app.py | 7 + pyvcell/vcml/models_math.py | 34 ++ pyvcell/vcml/vcml_reader.py | 38 ++- pyvcell/vcml/vcml_writer.py | 46 ++- tests/conftest.py | 1 + .../SmallSpatialProject_3D_well_mixed.vcml | 311 ++++++++++++++++++ tests/fixtures/vcell_model_fixtures.py | 7 + tests/sim_results/test_region_variables.py | 71 ++++ tests/vcml/test_math_description.py | 33 +- 11 files changed, 582 insertions(+), 4 deletions(-) create mode 100644 tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml create mode 100644 tests/sim_results/test_region_variables.py diff --git a/pyvcell/_internal/simdata/zarr_writer.py b/pyvcell/_internal/simdata/zarr_writer.py index 090cc4d..3f8bc93 100644 --- a/pyvcell/_internal/simdata/zarr_writer.py +++ b/pyvcell/_internal/simdata/zarr_writer.py @@ -20,10 +20,20 @@ def write_zarr(pde_dataset: PdeDataSet, data_functions: DataFunctions, mesh: Car volume_data_vars: list[DataBlockHeader] = [ v for v in pde_dataset.variables_block_headers() if v.var_info.variable_type == VariableType.VOLUME ] + # Volume-region state variables (well-mixed species: one value per connected volume region, as + # ``::``) are replicated onto every voxel of their region, so they read like volume + # variables — a channel each, and a numexpr binding for the functions that use them. The solver's + # built-in region quantities (``vcRegionVolume``, no domain) are not state and stay out. + volume_region_vars: list[DataBlockHeader] = [ + v + for v in pde_dataset.variables_block_headers() + if v.var_info.variable_type == VariableType.VOLUME_REGION and "::" in v.var_info.var_name + ] volume_functions: list[NamedFunction] = [ f for f in data_functions.named_functions if f.variable_type == VariableType.VOLUME ] - num_channels = len(volume_data_vars) + len(volume_functions) + 5 # 5 extra channels for region map, t, x, y, z + # 5 extra channels for region map, t, x, y, z + num_channels = len(volume_data_vars) + len(volume_region_vars) + len(volume_functions) + 5 num_t: int = len(pde_dataset.times()) times: list[float] = pde_dataset.times() header = pde_dataset.first_data_zip_file_metadata().file_header @@ -142,6 +152,28 @@ def _get_domain_mask(dname: str) -> np.ndarray: channel_metadata[c]["mean_values"].append(float(np.mean(masked))) c = c + 1 + # add volume-region state variables, replicated per region onto the voxels + for v in volume_region_vars: + region_values: NDArray[np.float64] = np.asarray(pde_dataset.get_data(v.var_info, times[t])) + var_data = region_values[region_map_flat].reshape((num_z, num_y, num_x)) + z1[t, c, :, :, :] = var_data + domain_name, var_name = v.var_info.var_name.split("::") + bindings[var_name] = var_data + if t == 0: + channel_metadata.append({ + "index": c, + "label": var_name, + "domain_name": domain_name, + "min_values": [], + "max_values": [], + "mean_values": [], + }) + masked = var_data[_get_domain_mask(domain_name)] + channel_metadata[c]["min_values"].append(float(np.min(masked))) + channel_metadata[c]["max_values"].append(float(np.max(masked))) + channel_metadata[c]["mean_values"].append(float(np.mean(masked))) + c = c + 1 + # add volumetric functions for f in volume_functions: evaluated_data = f.evaluate(variable_bindings=bindings) diff --git a/pyvcell/vcml/__init__.py b/pyvcell/vcml/__init__.py index dda84f9..c4b9148 100644 --- a/pyvcell/vcml/__init__.py +++ b/pyvcell/vcml/__init__.py @@ -56,6 +56,7 @@ MathFunction, MathVariable, MathVariableType, + MembraneRegionEquation, MembraneSubDomain, OdeEquation, ParticleInitialCount, @@ -64,6 +65,7 @@ PdeEquation, VariableInitialCount, Velocity, + VolumeRegionEquation, ) from pyvcell.vcml.vcml_reader import VcmlReader @@ -211,6 +213,7 @@ def __dir__() -> list[str]: "MathFunction", "MathVariable", "MathVariableType", + "MembraneRegionEquation", "MembraneSubDomain", "Model", "ModelParameter", @@ -240,6 +243,7 @@ def __dir__() -> list[str]: "VcmlWriter", "Velocity", "Version", + "VolumeRegionEquation", "cartesian_mesh_from_geometry", "connect", "field_data_refs", diff --git a/pyvcell/vcml/models_app.py b/pyvcell/vcml/models_app.py index 41953b8..1802d30 100644 --- a/pyvcell/vcml/models_app.py +++ b/pyvcell/vcml/models_app.py @@ -97,6 +97,13 @@ class SpeciesMapping(VcmlNode): velocity_y: float | str | None = None velocity_z: float | str | None = None boundary_values: list[float | str | None] = Field(default_factory=list) + # VCell's LocalizedCompoundSpec flags. ``force_constant``: a clamped species (its initial condition for + # all time). ``well_mixed``: in a spatial application, one value per region — the math becomes a + # VolumeRegionVariable / MembraneRegionVariable. ``force_continuous``: keep a species deterministic in a + # hybrid stochastic application. + force_constant: bool = False + well_mixed: bool = False + force_continuous: bool = False @property def expressions(self) -> list[str]: diff --git a/pyvcell/vcml/models_math.py b/pyvcell/vcml/models_math.py index 6e2badd..34664c3 100644 --- a/pyvcell/vcml/models_math.py +++ b/pyvcell/vcml/models_math.py @@ -108,6 +108,38 @@ class PdeEquation(VcmlNode): velocity: Velocity | None = None +class VolumeRegionEquation(VcmlNode): + """An ODE for a volume-region variable: one value per connected region of the compartment. + + ``d/dt u_R = uniform_rate + (1/|R|) ∫_R volume_rate dV`` over each region ``R`` (plus the fluxes the + membranes bounding ``R`` contribute through their jump conditions). ``solution_type`` is VCell's + ```` ("unknown" or "exact"); an exact solution is ``solution``. + """ + + name: str + uniform_rate: str | None = None + volume_rate: str | None = None + initial: str | None = None + solution_type: str | None = None + solution: str | None = None + + +class MembraneRegionEquation(VcmlNode): + """An ODE for a membrane-region variable: one value per connected region of the membrane. + + ``d/dt u_R = uniform_rate + (1/|R|) ∫_R membrane_rate dA`` over each membrane region ``R`` — most often + the membrane potential of an electrically coupled model. ``solution_type`` / ``solution`` as for + :class:`VolumeRegionEquation`. + """ + + name: str + uniform_rate: str | None = None + membrane_rate: str | None = None + initial: str | None = None + solution_type: str | None = None + solution: str | None = None + + class JumpCondition(VcmlNode): """A membrane flux jump condition for a volume variable.""" @@ -181,6 +213,7 @@ class CompartmentSubDomain(VcmlNode): boundary_types: list[MathBoundaryType] = Field(default_factory=list) ode_equations: list[OdeEquation] = Field(default_factory=list) pde_equations: list[PdeEquation] = Field(default_factory=list) + volume_region_equations: list[VolumeRegionEquation] = Field(default_factory=list) variable_initial_counts: list[VariableInitialCount] = Field(default_factory=list) jump_processes: list[JumpProcess] = Field(default_factory=list) particle_jump_processes: list[ParticleJumpProcess] = Field(default_factory=list) @@ -196,6 +229,7 @@ class MembraneSubDomain(VcmlNode): boundary_types: list[MathBoundaryType] = Field(default_factory=list) ode_equations: list[OdeEquation] = Field(default_factory=list) pde_equations: list[PdeEquation] = Field(default_factory=list) + membrane_region_equations: list[MembraneRegionEquation] = Field(default_factory=list) jump_conditions: list[JumpCondition] = Field(default_factory=list) particle_jump_processes: list[ParticleJumpProcess] = Field(default_factory=list) particle_properties: list[ParticleProperties] = Field(default_factory=list) diff --git a/pyvcell/vcml/vcml_reader.py b/pyvcell/vcml/vcml_reader.py index 3d24d29..549badd 100644 --- a/pyvcell/vcml/vcml_reader.py +++ b/pyvcell/vcml/vcml_reader.py @@ -339,6 +339,12 @@ def _parse_compartment_subdomain(self, element: _Element) -> vcm.CompartmentSubD subdomain.ode_equations.append(self._parse_ode_equation(child)) elif tag == "PdeEquation": subdomain.pde_equations.append(self._parse_pde_equation(child)) + elif tag == "VolumeRegionEquation": + subdomain.volume_region_equations.append( + vcm.VolumeRegionEquation.model_validate( + self._parse_region_equation(child, "VolumeRate", "volume_rate") + ) + ) elif tag in ("VariableInitialCount", "VariableInitialPoissonExpectedCount"): subdomain.variable_initial_counts.append( vcm.VariableInitialCount( @@ -368,6 +374,12 @@ def _parse_membrane_subdomain(self, element: _Element) -> vcm.MembraneSubDomain: subdomain.ode_equations.append(self._parse_ode_equation(child)) elif tag == "PdeEquation": subdomain.pde_equations.append(self._parse_pde_equation(child)) + elif tag == "MembraneRegionEquation": + subdomain.membrane_region_equations.append( + vcm.MembraneRegionEquation.model_validate( + self._parse_region_equation(child, "MembraneRate", "membrane_rate") + ) + ) elif tag == "JumpCondition": subdomain.jump_conditions.append(self._parse_jump_condition(child)) elif tag in ("ParticleJumpProcess", "LangevinParticleJumpProcess"): @@ -376,6 +388,25 @@ def _parse_membrane_subdomain(self, element: _Element) -> vcm.MembraneSubDomain: subdomain.particle_properties.append(self._parse_particle_properties(child)) return subdomain + def _parse_region_equation(self, element: _Element, rate_tag: str, rate_field: str) -> dict[str, str | None]: + """The fields of a ``VolumeRegionEquation`` / ``MembraneRegionEquation``: ``UniformRate``, the + region's own rate (``rate_tag``), ``Initial``, and ```` whose + text is the exact solution.""" + fields: dict[str, str | None] = {"name": element.get("Name", default="unnamed")} + for child in element: + tag = strip_namespace(child.tag) + if tag == "UniformRate": + fields["uniform_rate"] = _text(child) + elif tag == rate_tag: + fields[rate_field] = _text(child) + elif tag == "Initial": + fields["initial"] = _text(child) + elif tag == "SolutionType": + fields["solution_type"] = child.get("Type") + if child.get("Type") == "exact": + fields["solution"] = _text(child) + return fields + def _parse_ode_equation(self, element: _Element) -> vcm.OdeEquation: equation = vcm.OdeEquation( name=element.get("Name", default="unnamed"), solution_type=element.get("SolutionType") @@ -613,7 +644,12 @@ def visit_BoundariesTypes(self, element: _Element, node: vc.CompartmentMapping) def visit_LocalizedCompoundSpec(self, element: _Element, node: vc.Application) -> None: species_name: str = element.get("LocalizedCompoundRef", default="unnamed") - species_mapping = vc.SpeciesMapping(species_name=species_name) + species_mapping = vc.SpeciesMapping( + species_name=species_name, + force_constant=element.get("ForceConstant") == "true", + well_mixed=element.get("WellMixed") == "true", + force_continuous=element.get("ForceContinuous") == "true", + ) node.species_mappings.append(species_mapping) self.generic_visit_children(element, species_mapping) diff --git a/pyvcell/vcml/vcml_writer.py b/pyvcell/vcml/vcml_writer.py index 4423620..b1f06f1 100644 --- a/pyvcell/vcml/vcml_writer.py +++ b/pyvcell/vcml/vcml_writer.py @@ -26,11 +26,13 @@ JumpProcess, MathBoundaryType, MathDescription, + MembraneRegionEquation, MembraneSubDomain, OdeEquation, ParticleJumpProcess, ParticleProperties, PdeEquation, + VolumeRegionEquation, ) @@ -218,7 +220,13 @@ def write_application(self, application: Application, parent: _Element) -> None: reaction_context_element = Element("ReactionContext") parent.append(reaction_context_element) for species_mapping in application.species_mappings: - mapping_element = Element("LocalizedCompoundSpec", LocalizedCompoundRef=species_mapping.species_name) + mapping_element = Element( + "LocalizedCompoundSpec", + LocalizedCompoundRef=species_mapping.species_name, + ForceConstant=str(species_mapping.force_constant).lower(), + WellMixed=str(species_mapping.well_mixed).lower(), + ForceContinuous=str(species_mapping.force_continuous).lower(), + ) reaction_context_element.append(mapping_element) self.write_species_mapping(species_mapping, mapping_element) for reaction_mapping in application.reaction_mappings: @@ -466,6 +474,14 @@ def write_compartment_subdomain(self, subdomain: CompartmentSubDomain, parent: _ self.write_ode_equation(ode_equation, subdomain_element) for pde_equation in subdomain.pde_equations: self.write_pde_equation(pde_equation, subdomain_element) + for volume_region_equation in subdomain.volume_region_equations: + self.write_region_equation( + "VolumeRegionEquation", + volume_region_equation, + "VolumeRate", + volume_region_equation.volume_rate, + subdomain_element, + ) for variable_initial_count in subdomain.variable_initial_counts: tag = "VariableInitialPoissonExpectedCount" if variable_initial_count.poisson else "VariableInitialCount" count_element = Element(tag, Name=variable_initial_count.name) @@ -491,6 +507,14 @@ def write_membrane_subdomain(self, subdomain: MembraneSubDomain, parent: _Elemen self.write_ode_equation(ode_equation, subdomain_element) for pde_equation in subdomain.pde_equations: self.write_pde_equation(pde_equation, subdomain_element) + for membrane_region_equation in subdomain.membrane_region_equations: + self.write_region_equation( + "MembraneRegionEquation", + membrane_region_equation, + "MembraneRate", + membrane_region_equation.membrane_rate, + subdomain_element, + ) for jump_condition in subdomain.jump_conditions: self.write_jump_condition(jump_condition, subdomain_element) for particle_jump_process in subdomain.particle_jump_processes: @@ -542,6 +566,26 @@ def write_pde_equation(self, equation: PdeEquation, parent: _Element) -> None: equation_element.append(Element("Velocity", attrib=velocity_attrs)) self._append_text_element(equation_element, "Solution", equation.solution) + def write_region_equation( + self, + tag: str, + equation: VolumeRegionEquation | MembraneRegionEquation, + rate_tag: str, + rate: str | None, + parent: _Element, + ) -> None: + """A ``VolumeRegionEquation`` / ``MembraneRegionEquation`` as VCell's Xmlproducer writes it: the + rates default to ``0.0``; ```` carries an exact solution as its text.""" + equation_element = Element(tag, Name=equation.name) + parent.append(equation_element) + self._append_text_element(equation_element, "UniformRate", equation.uniform_rate or "0.0") + self._append_text_element(equation_element, rate_tag, rate or "0.0") + self._append_text_element(equation_element, "Initial", equation.initial) + solution_type = Element("SolutionType", Type=equation.solution_type or "unknown") + if solution_type.get("Type") == "exact" and equation.solution is not None: + solution_type.text = equation.solution + equation_element.append(solution_type) + def write_jump_condition(self, condition: JumpCondition, parent: _Element) -> None: condition_element = Element("JumpCondition", Name=condition.name) parent.append(condition_element) diff --git a/tests/conftest.py b/tests/conftest.py index 4ba4c4e..38587f5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -29,6 +29,7 @@ vcml_spatial_model_1d_path, vcml_spatial_particle_path, vcml_spatial_small_3d_path, + vcml_spatial_well_mixed_3d_path, vcml_tutorial_multiapp_pde_path, ) diff --git a/tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml b/tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml new file mode 100644 index 0000000..554a68a --- /dev/null +++ b/tests/fixtures/data/SmallSpatialProject_3D_well_mixed.vcml @@ -0,0 +1,311 @@ + + + + + + + 1.0 + 0.5 + 1.0 + 1.0 + 1.0 + 1.0 + + + s0 + + + s1 + + + s2 + + + s3 + + + + + + + + + + + + + ((Kf_r0 * s0) - (1.0 * Kr_r0 * s1)) + + + + + + + ((Kf_r1 * s0) - (Kr_r1 * s2)) + 0.0 + 1.0 + + + + + + + ((Kf_r2 * s2) - (Kr_r2 * s3)) + 0.0 + 1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ((pow(( - 5.0 + x),2.0) + pow(( - 5.0 + y),2.0) + pow(( - 5.0 + z),2.0)) < 16.0) + + + 1.0 + + + + + + + + + + + + + + + + + + + + + + (1.0 + sin(x)) + 1.0E-4 + + + + (1.0 + cos(x)) + 1.0E-4 + + + + (1.0 + (sin(x) * cos(y))) + 1.0E-4 + + + + (1.0 + cos(y)) + 1.0000000000000002E-6 + + + + + + + + 96485.3321 + 9.64853321E-5 + 1.0E-9 + 6.02214179E11 + 3.141592653589793 + 8314.46261815 + 300.0 + 1.0 + 1000.0 + 1.0 + 1.0 + 1.0 + 0.001660538783162726 + 0.5 + 1.0 + 1.0 + 1.0 + 1.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 1.0E-4 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 0.0 + 1.0E-4 + 0.0 + 1.0 + 1.0 + + + + + ((Kf_r0 * s0) - (1.0 * Kr_r0 * s1)) + ((Kf_r1 * s0) - (Kr_r1 * s2)) + ((Kf_r2 * s2) - (Kr_r2 * s3)) + (AreaPerUnitArea_m0 / VolumePerUnitVolume_c0) + (AreaPerUnitArea_m0 / VolumePerUnitVolume_c1) + (1.0 + sin(x)) + (1.0 + cos(x)) + (1.0 + cos(y)) + (1.0 + (sin(x) * cos(y))) + (VolumePerUnitVolume_c0 * vcRegionVolume('subdomain0')) + (VolumePerUnitVolume_c1 * vcRegionVolume('subdomain1')) + (AreaPerUnitArea_m0 * vcRegionArea('subdomain0_subdomain1_membrane')) + vcRegionArea('subdomain0_subdomain1_membrane') + vcRegionVolume('subdomain0') + vcRegionVolume('subdomain1') + + + + + + + + + + 0.0 + s3_diffusionRate + s3_init_umol_l_1 + + + + + + + + + + + + - J_r0 + s0_diffusionRate + s0_init_umol_l_1 + + + 0.0 + J_r0 + s1_init_umol_l_1 + + + + + + + + + + + + 0.0 + (J_r1 - J_r2) + s2_init_umol_dm_2 + + + + ( - 1.0 * KFlux_m0_c0 * J_r1) + 0.0 + + + 0.0 + 0.0 + + + 0.0 + (KFlux_m0_c1 * J_r2) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/fixtures/vcell_model_fixtures.py b/tests/fixtures/vcell_model_fixtures.py index d6e063c..526fe29 100644 --- a/tests/fixtures/vcell_model_fixtures.py +++ b/tests/fixtures/vcell_model_fixtures.py @@ -19,6 +19,13 @@ def vcml_spatial_small_3d_path() -> Path: return FIXTURE_DATA_DIR / "SmallSpatialProject_3D.vcml" +@pytest.fixture +def vcml_spatial_well_mixed_3d_path() -> Path: + """SmallSpatialProject_3D with s1 (volume) and s2 (membrane) well-mixed, math regenerated by + libvcell: s1 becomes a VolumeRegionVariable, s2 a MembraneRegionVariable.""" + return FIXTURE_DATA_DIR / "SmallSpatialProject_3D_well_mixed.vcml" + + @pytest.fixture def vcml_spatial_bunny_3d_path() -> Path: return FIXTURE_DATA_DIR / "Bunny.vcml" diff --git a/tests/sim_results/test_region_variables.py b/tests/sim_results/test_region_variables.py new file mode 100644 index 0000000..b7d3973 --- /dev/null +++ b/tests/sim_results/test_region_variables.py @@ -0,0 +1,71 @@ +"""Well-mixed species in a spatial application: region variables through math, VCML and results. + +A well-mixed species has one value per connected region. VCell generates it as a VolumeRegionVariable / +MembraneRegionVariable with a region equation, and fvsolver writes it as one value per region. The zarr +writer replicates a volume-region variable onto every voxel of its region, so it reads like a volume +variable: a channel of its own, and a binding for the functions that use it. +""" + +from pathlib import Path + +import numpy as np + +import pyvcell.vcml as vc +from pyvcell._internal.simdata.simdata_models import VariableType + + +def test_species_mapping_flags_roundtrip(vcml_spatial_well_mixed_3d_path: Path) -> None: + """WellMixed / ForceConstant / ForceContinuous survive a read -> write -> read (they used to be dropped, + so a well-mixed species came back spatial, and a clamped species came back unclamped).""" + biomodel = vc.VcmlReader.biomodel_from_file(vcml_spatial_well_mixed_3d_path) + flags = {m.species_name: m.well_mixed for m in biomodel.applications[0].species_mappings} + assert flags == {"s0": False, "s1": True, "s2": True, "s3": False} + + mapping = next(m for m in biomodel.applications[0].species_mappings if m.species_name == "s0") + mapping.force_constant = True + vcml_str = vc.VcmlWriter().write_vcml(document=vc.VCMLDocument(biomodel=biomodel)) + assert 'LocalizedCompoundRef="s1" ForceConstant="false" WellMixed="true"' in vcml_str + roundtripped = vc.VcmlReader.biomodel_from_str(vcml_str) + assert roundtripped == biomodel + + +def test_volume_region_variable_in_results(vcml_spatial_well_mixed_3d_path: Path) -> None: + biomodel = vc.load_vcml_file(vcml_spatial_well_mixed_3d_path) + biomodel.applications[0].add_sim(name="sim", duration=1.0, output_time_step=0.25, mesh_size=(12, 12, 12)) + result = vc.simulate(biomodel, "sim") + + labels = [channel.label for channel in result.channel_data] + # s1 (well-mixed, subdomain0) is a channel after the volume variables; J_r0 = f(s0, s1) now evaluates + assert labels == [ + "region_mask", + "t", + "x", + "y", + "z", + "s0", + "s3", + "s1", + "J_r0", + "s0_init_uM", + "s1_init_uM", + "s3_init_uM", + ] + + header = next(h for h in result.pde_dataset.variables_block_headers() if h.var_info.var_name == "subdomain0::s1") + assert header.var_info.variable_type == VariableType.VOLUME_REGION + + zarr_data = result.zarr_dataset + region_map = np.asarray(zarr_data[0, 0]) + in_subdomain0 = np.isin(region_map, list(result.mesh.get_volume_region_ids("subdomain0"))) + s1_index = labels.index("s1") + times = result.pde_dataset.times() + for t_index, t in enumerate(times): + s1 = np.asarray(zarr_data[t_index, s1_index]) + per_region = np.asarray(result.pde_dataset.get_data(header.var_info, t)) + # every voxel holds its region's value: constant across subdomain0, equal to the solver's value + assert np.allclose(s1, per_region[region_map.astype(int)]) + assert np.ptp(s1[in_subdomain0]) == 0.0 + # the species actually evolves (so the replication is not trivially constant in time) + assert not np.isclose( + np.asarray(zarr_data[0, s1_index])[in_subdomain0][0], np.asarray(zarr_data[-1, s1_index])[in_subdomain0][0] + ) diff --git a/tests/vcml/test_math_description.py b/tests/vcml/test_math_description.py index 402f848..2392082 100644 --- a/tests/vcml/test_math_description.py +++ b/tests/vcml/test_math_description.py @@ -57,6 +57,35 @@ def test_math_description_membrane_jump_condition(vcml_spatial_bunny_3d_path: Pa assert jump_condition.out_flux is not None +def test_math_description_region_equations(vcml_spatial_well_mixed_3d_path: Path) -> None: + """Well-mixed species in a spatial application become region variables with region equations.""" + biomodel = vc.VcmlReader.biomodel_from_file(vcml_spatial_well_mixed_3d_path) + math = biomodel.applications[0].math_description + assert math is not None + var_types = {v.name: v.var_type for v in math.variables} + assert var_types["s1"] == vc.MathVariableType.volume_region + assert var_types["s2"] == vc.MathVariableType.membrane_region + + (subdomain0,) = (c for c in math.compartment_subdomains if c.name == "subdomain0") + assert subdomain0.volume_region_equations == [ + vc.VolumeRegionEquation( + name="s1", uniform_rate="0.0", volume_rate="J_r0", initial="s1_init_umol_l_1", solution_type="unknown" + ) + ] + assert "s1" not in [eq.name for eq in subdomain0.pde_equations + subdomain0.ode_equations] + + (membrane,) = math.membrane_subdomains + assert membrane.membrane_region_equations == [ + vc.MembraneRegionEquation( + name="s2", + uniform_rate="0.0", + membrane_rate="(J_r1 - J_r2)", + initial="s2_init_umol_dm_2", + solution_type="unknown", + ) + ] + + def test_math_description_ode_nonspatial(vcml_nonspatial_ode_path: Path) -> None: biomodel = vc.VcmlReader.biomodel_from_file(vcml_nonspatial_ode_path) math = biomodel.applications[0].math_description @@ -115,10 +144,11 @@ def test_math_description_roundtrip( vcml_spatial_particle_path: Path, vcml_spatial_model_1d_path: Path, vcml_spatial_bunny_3d_path: Path, + vcml_spatial_well_mixed_3d_path: Path, ) -> None: """Full biomodel write -> read preserves the math description (and everything else). - Covers nonspatial (ODE / stochastic) and spatial (PDE / membrane / particle) + Covers nonspatial (ODE / stochastic) and spatial (PDE / membrane / particle / region) geometries. """ paths = [ @@ -127,6 +157,7 @@ def test_math_description_roundtrip( vcml_spatial_particle_path, vcml_spatial_model_1d_path, vcml_spatial_bunny_3d_path, + vcml_spatial_well_mixed_3d_path, ] for path in paths: biomodel = vc.VcmlReader.biomodel_from_file(path)