diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index b9ff80d5..9541779f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -11,7 +11,8 @@ repos:
exclude: |
(?x)^(
joss-paper/JOSS_Fig2.png|
- examples/example2-3d/spine_mesh.xml
+ examples/example2-3d/spine_mesh.xml|
+ examples/example1/spheroid_ellipsoid_mesh.h5
)$
- id: check-docstring-first
- id: debug-statements
diff --git a/examples/example1/example1_CH.ipynb b/examples/example1/example1_CH.ipynb
new file mode 100644
index 00000000..08f60a2a
--- /dev/null
+++ b/examples/example1/example1_CH.ipynb
@@ -0,0 +1,685 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "f65f18d7",
+ "metadata": {},
+ "source": [
+ "# Example 1-CH: Cahn Hilliard patterns for 2D aggregation-diffusion\n",
+ "\n",
+ "In this case, we consider a simple 2D geometry comprised of two compartments:\n",
+ "- surf - 2D surface\n",
+ "- edge - outer edges of the surface (1D)\n",
+ "\n",
+ "We implement a Cahn-Hilliard model with two species, one exhibiting aggregation-diffusion ($B$) and the other purely diffusive ($X$).\n",
+ "The equations governing their evolution are given by:\n",
+ "\n",
+ "$$\n",
+ "\\partial_t{u_X} = -k_{on} u_X + k_{off} u_B + D_X \\nabla \\cdot (\\hat{\\mu}_X \\nabla u_X) \\\\\n",
+ "\\partial_t{u_B} = k_{on} u_X - k_{off} u_B + D_B \\nabla \\cdot (\\hat{\\mu}_B \\nabla u_B),\n",
+ "$$\n",
+ "\n",
+ "where aggregation is accounted for through the chemical potential of $B$, whose nondimensional version ($\\hat{\\mu}_B = \\frac{\\mu_B}{k_B T}$) is given by\n",
+ "\n",
+ "$$\n",
+ "\\hat{\\mu}_B = (\\ln \\phi_B - \\ln (1-\\phi_B)) - \\hat{A} (2 \\phi_B - 1) - \\frac{\\hat{A}}{u_{B,max}} \\nabla^2 \\phi_B\n",
+ "$$\n",
+ "\n",
+ "in which $\\phi_B$ is the area fraction occupied by species B, which is proportional to its concentration; that is, $\\phi_B = \\frac{u_B}{u_{B,max}}$.\n",
+ "The chemical potential is defined similarly for X, but with $\\hat{A}=0$; that is, $\\hat{\\mu}_X = (\\ln \\phi_X - \\ln (1-\\phi_X))$\n",
+ "\n",
+ "We solve these equations over a square domain with no-flux boundary conditions."
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "956a0fd1",
+ "metadata": {},
+ "source": [
+ "We begin with the necessary imports:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc398816",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import dolfin as d\n",
+ "import sympy as sym\n",
+ "import numpy as np\n",
+ "import pathlib\n",
+ "import gmsh # must be imported before pyvista if dolfin is imported first\n",
+ "\n",
+ "from smart import config, common, mesh, model, mesh_tools, visualization\n",
+ "from smart.units import unit\n",
+ "from smart.model_assembly import (\n",
+ " Compartment,\n",
+ " Parameter,\n",
+ " Reaction,\n",
+ " Species,\n",
+ " SpeciesContainer,\n",
+ " ParameterContainer,\n",
+ " CompartmentContainer,\n",
+ " ReactionContainer,\n",
+ ")\n",
+ "import logging\n",
+ "from matplotlib import pyplot as plt"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "028bb85e",
+ "metadata": {},
+ "source": [
+ "We will set the logging level to `INFO`. This will display some output during the simulation. If you want to get even more output you could set the logging level to `DEBUG`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c6e826d7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "logger = logging.getLogger(\"smart\")\n",
+ "logger.setLevel(logging.INFO)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "c7fd3be3",
+ "metadata": {},
+ "source": [
+ "Futhermore, you could also save the logs to a file by attaching a file handler to the logger as follows.\n",
+ "\n",
+ "```\n",
+ "file_handler = logging.FileHandler(\"filename.log\")\n",
+ "file_handler.setFormatter(logging.Formatter(smart.config.base_format))\n",
+ "logger.addHandler(file_handler)\n",
+ "```"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "95b9d865",
+ "metadata": {},
+ "source": [
+ "We define the various units for use in the model. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f4023cf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Aliases - base units\n",
+ "um = unit.um\n",
+ "molecule = unit.molecule\n",
+ "sec = unit.sec\n",
+ "dimensionless = unit.dimensionless\n",
+ "D_unit = um**2 / sec\n",
+ "flux_unit = molecule / (um * sec)\n",
+ "surf_unit = molecule / um**2\n",
+ "edge_unit = molecule / um"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "46582d26",
+ "metadata": {},
+ "source": [
+ "## Generate model\n",
+ "\n",
+ "### Compartments\n",
+ "As described above, the two compartments are the \"surf\" (2D) and edge (1D). These are initialized by calling:\n",
+ "```\n",
+ "compartment_var = Compartment(name, dimensionality, compartment_units, cell_marker)\n",
+ "```\n",
+ "where\n",
+ "- name: string naming the compartment\n",
+ "- dimensionality: topological dimensionality (e.g. 2 for surf, 1 for edge)\n",
+ "- compartment_units: length units for the compartment (um for both here)\n",
+ "- cell_marker: integer marker value identifying each compartment in the parent mesh"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "09079b17",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "surf = Compartment(\"surf\", 2, um, 10)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "2db8daf9",
+ "metadata": {},
+ "source": [
+ "Now we initialize a compartment container and add both compartments to it."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc3393cb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cc = CompartmentContainer()\n",
+ "cc.add([surf])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "8ee2332b",
+ "metadata": {},
+ "source": [
+ "### Species\n",
+ "In this case, we have a two species, \"X\" and \"B\", which exist in the 2D \"surf\" domain. Each is initialized by calling:\n",
+ "```\n",
+ "species_var = Species(\n",
+ " name, initial_condition, concentration_units,\n",
+ " D, diffusion_units, compartment_name, group (opt)\n",
+ " )\n",
+ "```\n",
+ "where\n",
+ "- name: string naming the species\n",
+ "- initial_condition: initial concentration for this species (can be an expression given by a string to be parsed by sympy - the only unknowns in the expression should be x, y, and z)\n",
+ "- concentration_units: concentration units for this species (molecules/μm2 here)\n",
+ "- D: diffusion coefficient\n",
+ "- diffusion_units: units for diffusion coefficient (μm2/sec here)\n",
+ "- compartment_name: each species should be assigned to a single compartment (\"surf\", here)\n",
+ "- group (opt): for larger models, specifies a group of species this belongs to;\n",
+ " for organizational purposes when there are multiple reaction modules\n",
+ "\n",
+ "With the added CH features, we also must provide the following for Cahn-Hilliard type species:\n",
+ "- `CH = True` - this tells SMART that we are considering both aggregation and diffusion\n",
+ "- umax: maximum surface density of X or B\n",
+ "- A_hat: strength of aggregation\n",
+ "\n",
+ "Note that A_hat is dimensionless here and the chemical potential is a variable generated *internally* within this branch of SMART. Because only the nondimensional chemical potential appears in the dynamical equation for $u_A$, the chemical potential is always normalized to the thermal energy scale $k_B T$. (see equations up top for consistency on this point)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3f6f384b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "l0 = np.sqrt(4*np.pi) # reference length scale\n",
+ "Shat = 10.0\n",
+ "phi0_X = 0.1\n",
+ "phi0_B = 0.1\n",
+ "sigma_s = Shat/l0**2\n",
+ "X = Species(\"X\", phi0_X*sigma_s, surf_unit, 1.0, D_unit, \"surf\", CH=True, umax=sigma_s, A_hat = 0) \n",
+ "B = Species(\"B\", phi0_B*sigma_s, surf_unit, 1.0, D_unit, \"surf\", CH=True, umax=sigma_s, A_hat = 50)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "f77d2d31",
+ "metadata": {},
+ "source": [
+ "Create a species container and add both species to it:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5c1df887",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sc = SpeciesContainer()\n",
+ "sc.add([X, B])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "74d1d353",
+ "metadata": {},
+ "source": [
+ "### Parameters and Reactions\n",
+ "Parameters and reactions are generally defined together, although the order does not strictly matter. Parameters are specified as:\n",
+ "```\n",
+ "param_var = Parameter(name, value, unit, group (opt), notes (opt), use_preintegration (opt))\n",
+ "```\n",
+ "where\n",
+ "- name: string naming the parameter\n",
+ "- value: value of the given parameter\n",
+ "- unit: units associated with given value\n",
+ "- group (optional): optional string placing this reaction in a reaction group; for organizational purposes when there are multiple reaction modules\n",
+ "- notes (optional): string related to this parameter\n",
+ "- use_preintegration (optional): in the case of a time-dependent parameter, uses preintegration in the solution process\n",
+ "\n",
+ "Reactions are specified by a variable number of arguments (arguments are indicated by (opt) are either never\n",
+ "required or only required in some cases, for more details see notes below and API documentation):\n",
+ "```\n",
+ "reaction_var = Reaction(\n",
+ " name, lhs, rhs, param_map,\n",
+ " eqn_f_str (opt), eqn_r_str (opt), reaction_type (opt), species_map,\n",
+ " explicit_restriction_to_domain (opt), group (opt), flux_scaling (opt)\n",
+ " )\n",
+ "```\n",
+ "- name: string naming the reaction\n",
+ "- lhs: list of strings specifying the reactants for this reaction\n",
+ "- rhs: list of strings specifying the products for this reaction\n",
+ " ***NOTE: the lists \"reactants\" and \"products\" determine the stoichiometry of the reaction;\n",
+ " for instance, if two A's react to give one B, the reactants list would be [\"A\",\"A\"],\n",
+ " and the products list would be [\"B\"]\n",
+ "- param_map: relationship between the parameters specified in the reaction string and those given\n",
+ " in the parameter container. By default, the reaction parameters are \"kon\" and \"koff\" when\n",
+ " a system obeys simple mass action. If the forward rate is given by a parameter \"k1\" and the\n",
+ " reverse rate is given by \"k2\", then param_map = {\"on\":\"k1\", \"off\":\"k2\"}\n",
+ "- eqn_f_str: For systems not obeying simple mass action, this string specifies the forward reaction rate\n",
+ " By default, this string is \"on*{all reactants multiplied together}\"\n",
+ "- eqn_r_str: For systems not obeying simple mass action, this string specifies the reverse reaction rate\n",
+ " By default, this string is \"off*{all products multiplied together}\"\n",
+ "- reaction_type (opt): either \"custom\" or \"mass_action\" (default is \"mass_action\") [never a required argument]\n",
+ "- species_map: same format as param_map; required if other species not listed in reactants or products appear in the\n",
+ " reaction string\n",
+ "- explicit_restriction_to_domain: string specifying where the reaction occurs; required if the reaction is not\n",
+ " constrained by the reaction string (e.g., if production occurs only at the boundary,\n",
+ " as it does here, but the species being produced exists through the entire volume)\n",
+ "- group (opt): string placing this reaction in a reaction group; for organizational purposes when there are multiple reaction modules\n",
+ "- flux_scaling (opt): in certain cases, a given reactant or product may experience a scaled flux (for instance, if we assume that\n",
+ " some of the molecules are immediately sequestered after the reaction); in this case, to signify that this flux \n",
+ " should be rescaled, we specify ''flux_scaling = {scaled_species: scale_factor}'', where scaled_species is a\n",
+ " string specifying the species to be scaled and scale_factor is a number specifying the rescaling factor\n",
+ "\n",
+ "For this system, we do not define any reactions on the boundary (`edge`). This corresponds to assuming a no-flux boundary condition."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "df027853",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "kon_hat = 1.0\n",
+ "koff_hat = 1.0\n",
+ "tref = l0**2 / float(X.D)\n",
+ "kon = Parameter(\"kon\", kon_hat*sigma_s**(3/2)/tref, 1/sec)\n",
+ "koff = Parameter(\"koff\", koff_hat*sigma_s**(3/2)/tref, 1/sec)\n",
+ "# Conversion of X to B\n",
+ "r1 = Reaction(\"r1\", [\"X\"], [\"B\"],\n",
+ " param_map={\"kon\": \"kon\", \"koff\": \"koff\"},\n",
+ " eqn_f_str=\"X*kon - B*koff\")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "a0670e78",
+ "metadata": {},
+ "source": [
+ "Create parameter and reaction containers and add in associated objects."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1eb19fb6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pc = ParameterContainer()\n",
+ "pc.add([kon, koff])\n",
+ "rc = ReactionContainer()\n",
+ "rc.add([r1])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "15c35d39",
+ "metadata": {},
+ "source": [
+ "## Create/load in mesh\n",
+ "\n",
+ "In SMART we have different levels of meshes. Here we create a UnitSquare mesh defined by\n",
+ "\n",
+ "$$\n",
+ "\\Omega = [0, 1] \\times [0, 1] \\subset \\mathbb{R}^2\n",
+ "$$\n",
+ "\n",
+ "which will serve as our parent mesh\n",
+ "\n",
+ "For our two domains, we have two associated \"child meshes\", which are set by the marker functions `mf2` and `mf1`:\n",
+ "- surf: in this case, all cells (triangles) belong to this mesh; here, marked by `mf2 = 1`\n",
+ "- edge: 1D child mesh including all line elements along the edges of the domain; here, marked by `mf1 = 3`\n",
+ "\n",
+ "Note that the marker values must be chosen to match those given in the compartment definitions above."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fe56e162",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "useSpheroid = True\n",
+ "if useSpheroid:\n",
+ " # rOuter = [0.6849, 0.4365, 2.1896]\n",
+ " # rInner = [0.0,0.0,0.0]\n",
+ " # domain, facet_markers, cell_markers = mesh_tools.create_ellipsoids(rOuter, rInner, hEdge=0.05)\n",
+ " vol = Compartment(\"vol\", 3, um, 1) # SMART just needs to know its a 3d mesh\n",
+ " cc.add(vol)\n",
+ " mesh_file = pathlib.Path(\"spheroid_ellipsoid_mesh.h5\")\n",
+ " # mesh_file = pathlib.Path(\"spheroid_ellipsoid_mesh_new.h5\")\n",
+ " # mesh_tools.write_mesh(domain, facet_markers, cell_markers, filename=mesh_file)\n",
+ "else:\n",
+ " # define dimensions of domain\n",
+ " Shat = 200\n",
+ " x_size = np.sqrt(Shat/B.umax)\n",
+ " y_size = np.sqrt(Shat/B.umax)\n",
+ " # Create mesh\n",
+ " m = 30\n",
+ " n = int(x_size/y_size)*m\n",
+ " rect_mesh = d.RectangleMesh(d.Point(0.0, 0.0), d.Point(x_size, y_size), n, m)\n",
+ " mf2 = d.MeshFunction(\"size_t\", rect_mesh, 2, 10)\n",
+ " mf1 = d.MeshFunction(\"size_t\", rect_mesh, 1, 0)\n",
+ " class OuterEdge(d.SubDomain):\n",
+ " def inside(self, x, on_boundary):\n",
+ " return on_boundary\n",
+ " outerEdge = OuterEdge()\n",
+ " outerEdge.mark(mf1, 3)\n",
+ " mesh_folder = pathlib.Path(\"rect_mesh\")\n",
+ " mesh_folder.mkdir(exist_ok=True)\n",
+ " mesh_file = mesh_folder / \"rect_mesh.h5\"\n",
+ " mesh_tools.write_mesh(rect_mesh, mf1, mf2, mesh_file)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "173e5474",
+ "metadata": {},
+ "source": [
+ "Finally, we initialize the `mesh.ParentMesh` object, using the hdf5 file as input."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bcbd06cd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "parent_mesh = mesh.ParentMesh(\n",
+ " mesh_filename=str(mesh_file),\n",
+ " mesh_filetype=\"hdf5\",\n",
+ " name=\"parent_mesh\",\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "8f7ac819",
+ "metadata": {},
+ "source": [
+ "## Initialize model and solver\n",
+ "Now we are ready to set up the model. First we load the default configurations and set the solver config."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "43170d20",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "config_cur = config.Config()\n",
+ "config_cur.flags.update({\"allow_unused_components\": True})\n",
+ "config_cur.solver.update(\n",
+ " {\n",
+ " \"final_t\": 100.0,\n",
+ " \"initial_dt\": 0.001,\n",
+ " \"time_precision\": 8,\n",
+ " \"attempt_timestep_restart_on_divergence\": True,\n",
+ " }\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "96783316",
+ "metadata": {},
+ "source": [
+ "We create the model object initialize the model using the `initialize` function found in the `smart.model` module. We then save the model information to a .pkl file for later reference.\n",
+ "\n",
+ "Note that we could later load the model information from the pickle file using the line:\n",
+ "```\n",
+ "model_cur = model.from_pickle(model_cur.pkl)\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a86c7435",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n",
+ "model_cur.initialize()\n",
+ "model_cur.to_pickle('model_cur.pkl')"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "c0d464e3",
+ "metadata": {},
+ "source": [
+ "We then perturb the initial conditions by adding white noise to the dolfin vectors associated with each species."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4c42042c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# add white noise perturbation to initial conditions\n",
+ "for sp_str in (\"B\"):#(\"X\", \"B\"):\n",
+ " sp = model_cur.sc[sp_str]\n",
+ " u = model_cur.cc[sp.compartment_name].u[\"u\"]\n",
+ " indices = sp.dof_map\n",
+ " uvec = u.vector()\n",
+ " values = uvec.get_local()\n",
+ " cur_seed = ord(sp_str) # set seed for reproducibility\n",
+ " generator_cur = np.random.default_rng(cur_seed)\n",
+ " values[indices] = np.multiply(values[indices],\n",
+ " generator_cur.normal(1, 0.01, len(indices)))\n",
+ " uvec.set_local(values)\n",
+ " uvec.apply(\"insert\")\n",
+ " nvec = model_cur.cc[sp.compartment_name].u[\"n\"].vector()\n",
+ " nvec.set_local(values)\n",
+ " nvec.apply(\"insert\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fbe592a7",
+ "metadata": {},
+ "source": [
+ "Define other functions to be used for mass conservation and cutoffs here"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8aa0f4cf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def _set_clipped_sum_c_tmp_scalar(sp, xi1_scalar, dt_val, epsilon):\n",
+ " \"\"\"clip(c + dt*xi_1) on all DOFs; xi_1 is a global scalar.\"\"\"\n",
+ " lo = sp.umax*(epsilon)\n",
+ " hi = sp.umax*(1.0 - epsilon)\n",
+ " c_tmp = d.Function(sp.V)\n",
+ " d.assign(c_tmp, sp.sol)\n",
+ " cvec = c_tmp.vector()[:]\n",
+ " cvec = np.clip(cvec + float(dt_val) * xi1_scalar, float(lo), float(hi))\n",
+ " c_tmp.vector().set_local(cvec)\n",
+ " c_tmp.vector().apply(\"insert\")\n",
+ " return c_tmp\n",
+ "\n",
+ "Xfunc = model_cur.sc[\"X\"].sol\n",
+ "Xdof = model_cur.sc[\"X\"].dof_map\n",
+ "Bfunc = model_cur.sc[\"B\"].sol\n",
+ "Bdof = model_cur.sc[\"B\"].dof_map\n",
+ "dx = d.Measure(\"dx\", model_cur.cc[\"surf\"].dolfin_mesh)\n",
+ "c_mass_init = d.assemble_mixed((Xfunc+Bfunc)*dx)\n",
+ "def F_mass(xi_arg1, epsilon, c_mass_init):\n",
+ " \"\"\"F(xi_1) = int clip(phi_X + dt*xi_1) + int clip(phi_B + dt*xi_1) minus previous X+B.\"\"\"\n",
+ " dt_val = float(model_cur.dt)\n",
+ " Xtemp = _set_clipped_sum_c_tmp_scalar(\n",
+ " model_cur.sc[\"X\"], xi_arg1, dt_val, epsilon)\n",
+ " Btemp = _set_clipped_sum_c_tmp_scalar(\n",
+ " model_cur.sc[\"B\"], xi_arg1, dt_val, epsilon)\n",
+ " dx = d.Measure(\"dx\", Xtemp.function_space().mesh())\n",
+ " mass_err = d.assemble_mixed((Xtemp+Btemp)*dx) - c_mass_init\n",
+ " return mass_err\n",
+ "\n",
+ "def assign_from_sub(subfunc, func, dofmap):\n",
+ " fullvec = func.vector()[:]\n",
+ " subvec = subfunc.vector()[:]\n",
+ " fullvec[dofmap] = subvec\n",
+ " func.vector().set_local(fullvec)\n",
+ " func.vector().apply(\"insert\")\n",
+ "\n",
+ "# check that the above functions are consistent\n",
+ "if F_mass(0.0, 0.01, c_mass_init) > 0.0:\n",
+ " raise ValueError(\"This should not be possible\")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "b610b5b8",
+ "metadata": {},
+ "source": [
+ "## Solve the system and write output data\n",
+ "Now, we are ready to start the solution process. We store the initial conditions to output files and then solve the system at each time step using the `monolithic_solve` function. Once we pass the final time chosen above, we exit the loop."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ce499ef3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Write initial condition(s) to file\n",
+ "results = dict()\n",
+ "result_folder = pathlib.Path(\"resultsRect\")\n",
+ "result_folder.mkdir(exist_ok=True)\n",
+ "for species_name, species in model_cur.sc.items:\n",
+ " results[species_name] = d.XDMFFile(\n",
+ " model_cur.mpi_comm_world, str(result_folder / f\"{species_name}.xdmf\")\n",
+ " )\n",
+ " results[species_name].parameters[\"flush_output\"] = True\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ "\n",
+ "# Set loglevel to warning in order not to pollute notebook output\n",
+ "logger.setLevel(logging.WARNING)\n",
+ "\n",
+ "epsilon = 0.001\n",
+ "xi_secant_max_iter = 500\n",
+ "\n",
+ "# Solve\n",
+ "while True:\n",
+ " print(f\"Time is {model_cur.t}\")\n",
+ " # Solve the system\n",
+ " model_cur.monolithic_solve()\n",
+ " model_cur.adjust_dt()\n",
+ " # Secant on global scalar xi_1: initial guesses (0, -dt).\n",
+ " xi_secant_iter = 0\n",
+ " xi_guess_prev = 0.0\n",
+ " xi_guess = -float(model_cur.dt)\n",
+ " secant_tol = 1e-12\n",
+ " F1 = F_mass(xi_guess, epsilon, c_mass_init)\n",
+ " F0 = F_mass(xi_guess_prev, epsilon, c_mass_init)\n",
+ " while (xi_secant_iter < xi_secant_max_iter and \n",
+ " abs(F1) > secant_tol and abs(F0) > secant_tol):\n",
+ " xi_secant_iter += 1\n",
+ " denom = F1 - F0\n",
+ " if abs(denom) < 1e-30 or xi_guess == xi_guess_prev:\n",
+ " break\n",
+ " xi_guess_next = xi_guess - F1 * (xi_guess - xi_guess_prev) / denom\n",
+ " xi_guess_prev = float(xi_guess)\n",
+ " xi_guess = float(xi_guess_next)\n",
+ " F1 = F_mass(xi_guess, epsilon, c_mass_init)\n",
+ " F0 = F_mass(xi_guess_prev, epsilon, c_mass_init)\n",
+ " print(f\"Secant approach converged in {xi_secant_iter} iterations\")\n",
+ " # now assign corrected values\n",
+ " Xnew = _set_clipped_sum_c_tmp_scalar(model_cur.sc[\"X\"], xi_guess, model_cur.dt, epsilon)\n",
+ " assign_from_sub(Xnew, Xfunc, Xdof)\n",
+ " Bnew = _set_clipped_sum_c_tmp_scalar(model_cur.sc[\"B\"], xi_guess, model_cur.dt, epsilon)\n",
+ " assign_from_sub(Bnew, Bfunc, Bdof)\n",
+ "\n",
+ "\n",
+ " for species_name, species in model_cur.sc.items:\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ " # End if we've passed the final time\n",
+ " if model_cur.t >= model_cur.final_t:\n",
+ " break"
+ ]
+ }
+ ],
+ "metadata": {
+ "jupytext": {
+ "cell_metadata_filter": "-all",
+ "main_language": "python",
+ "notebook_metadata_filter": "-all"
+ },
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/example1/spheroid_ellipsoid_mesh.h5 b/examples/example1/spheroid_ellipsoid_mesh.h5
new file mode 100644
index 00000000..41202e48
Binary files /dev/null and b/examples/example1/spheroid_ellipsoid_mesh.h5 differ
diff --git a/examples/example9/example9.ipynb b/examples/example9/example9.ipynb
new file mode 100644
index 00000000..4c6b7ccc
--- /dev/null
+++ b/examples/example9/example9.ipynb
@@ -0,0 +1,299 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "f65f18d7",
+ "metadata": {},
+ "source": [
+ "# Example 9: Advection-diffusion in the case of confined migration\n",
+ "Here, we consider the binding or uptake of a molecule in the case of confined migration. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc398816",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import dolfin as d\n",
+ "import sympy as sym\n",
+ "import numpy as np\n",
+ "import pathlib\n",
+ "import logging\n",
+ "import gmsh # must be imported before pyvista if dolfin is imported first\n",
+ "\n",
+ "from smart import config, mesh, model, mesh_tools, visualization\n",
+ "from smart.units import unit\n",
+ "from smart.model_assembly import (\n",
+ " Compartment,\n",
+ " Parameter,\n",
+ " Reaction,\n",
+ " Species,\n",
+ " SpeciesContainer,\n",
+ " ParameterContainer,\n",
+ " CompartmentContainer,\n",
+ " ReactionContainer,\n",
+ ")\n",
+ "\n",
+ "from matplotlib import pyplot as plt\n",
+ "import matplotlib.image as mpimg\n",
+ "from matplotlib import rcParams\n",
+ "\n",
+ "logger = logging.getLogger(\"smart\")\n",
+ "logger.setLevel(logging.INFO)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "95b9d865",
+ "metadata": {},
+ "source": [
+ "We define the relevant units here."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f4023cf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Aliases - base units\n",
+ "uM = unit.uM\n",
+ "um = unit.um\n",
+ "molecule = unit.molecule\n",
+ "sec = unit.sec\n",
+ "dimensionless = unit.dimensionless\n",
+ "# Aliases - units used in model\n",
+ "D_unit = um**2 / sec\n",
+ "flux_unit = uM * um / sec\n",
+ "vol_unit = uM\n",
+ "surf_unit = molecule / um**2"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "46582d26",
+ "metadata": {},
+ "source": [
+ "## Model generation\n",
+ "\n",
+ "We define the compartments and species first, with their respective containers."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "02a000f2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "EC = Compartment(\"EC\", 2, um, 1)\n",
+ "Cyto = Compartment(\"Cyto\", 2, um, 2)\n",
+ "Tube = Compartment(\"Tube\", 1, um, 10)\n",
+ "PM = Compartment(\"PM\", 1, um, 12)\n",
+ " # vel=[\"0\",\"0\",\"100.0*[1-(x[0]**2 + x[1]**2//4)]\"])\n",
+ "\n",
+ "cc = CompartmentContainer()\n",
+ "cc.add([EC, Cyto, Tube, PM])\n",
+ "\n",
+ "A = Species(\"A\", 10.0, vol_unit, 1.0, D_unit, \"EC\")\n",
+ "Abound = Species(\"Abound\", 0.1, surf_unit, 0.1, D_unit, \"PM\")\n",
+ "B = Species(\"B\", 10.0, vol_unit, 0.1, D_unit, \"EC\")\n",
+ "Bcyto = Species(\"Bcyto\", 1.0, vol_unit, 0.1, D_unit, \"Cyto\")\n",
+ "sc = SpeciesContainer()\n",
+ "sc.add([A, Abound, B, Bcyto])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "3c56e840",
+ "metadata": {},
+ "source": [
+ "Define parameters and reactions, then place in respective containers.\n",
+ "* r1: release of A from PM"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2e1f6882",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# reactions at PM\n",
+ "kon = Parameter(\"kon\", 0.0, flux_unit/vol_unit)\n",
+ "koff = Parameter(\"koff\", 0.0, 1/sec)\n",
+ "r1 = Reaction(\"r1\", [\"A\"], [\"Abound\"],\n",
+ " param_map={\"on\":\"kon\",\"off\":\"koff\"},\n",
+ " eqn_f_str=\"on*A - off*Abound\",\n",
+ " explicit_restriction_to_domain=\"PM\")\n",
+ "kin = Parameter(\"kin\", 0.01, flux_unit/vol_unit)\n",
+ "kout = Parameter(\"kout\", 0.01, flux_unit/vol_unit)\n",
+ "r2 = Reaction(\"r2\", [\"B\"], [\"Bcyto\"],\n",
+ " param_map={\"kin\":\"kin\",\"kout\":\"kout\"},\n",
+ " eqn_f_str=\"kin*B - kout*Bcyto\",\n",
+ " explicit_restriction_to_domain=\"PM\")\n",
+ "\n",
+ "pc = ParameterContainer()\n",
+ "pc.add([kon, koff, kin, kout])\n",
+ "rc = ReactionContainer()\n",
+ "rc.add([r1,r2])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "15c35d39",
+ "metadata": {},
+ "source": [
+ "## Create and load in mesh\n",
+ "\n",
+ "Here, we consider cells embedded in a cube mesh. The source cell is located at (0,0,0) and 8 other cells are spread equidistant through the mesh."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fe56e162",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tubeRad = 5.0\n",
+ "gapSize = 1.0\n",
+ "cellVol = 1200.0\n",
+ "hEdge = 0.4\n",
+ "hInnerEdge = 0.1\n",
+ "domain, facet_markers, cell_markers = mesh_tools.create_confined(tubeRad, gapSize, cellVol, hEdge, hInnerEdge)\n",
+ "# Write mesh and meshfunctions to file\n",
+ "mesh_folder = pathlib.Path(\"mesh\")\n",
+ "mesh_folder.mkdir(exist_ok=True)\n",
+ "mesh_path = mesh_folder / \"cyl_mesh.h5\"\n",
+ "mesh_tools.write_mesh(\n",
+ " domain, facet_markers, cell_markers, filename=mesh_path\n",
+ ")\n",
+ "parent_mesh = mesh.ParentMesh(\n",
+ " mesh_filename=str(mesh_path),\n",
+ " mesh_filetype=\"hdf5\",\n",
+ " name=\"parent_mesh\",\n",
+ ")\n",
+ "# visualization.plot_dolfin_mesh(domain, cell_markers, facet_markers)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "0943588e",
+ "metadata": {},
+ "source": [
+ "Initialize model and solver."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ac88bdec",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "config_cur = config.Config()\n",
+ "config_cur.flags.update({\"allow_unused_components\": True})\n",
+ "config_cur.flags.update({\"axisymmetric_model\": True})\n",
+ "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n",
+ "config_cur.solver.update(\n",
+ " {\n",
+ " \"final_t\": 1000.0,\n",
+ " \"initial_dt\": 0.01,\n",
+ " \"time_precision\": 8,\n",
+ " \"reset_timestep_for_negative_solution\": False,\n",
+ " }\n",
+ ")\n",
+ "model_cur.initialize()"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "5d5aacbd",
+ "metadata": {},
+ "source": [
+ "Initialize XDMF files for saving results, save model information to .pkl file, then solve the system until `model_cur.t > model_cur.final_t`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b54d28ca",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Write initial condition(s) to file\n",
+ "results = dict()\n",
+ "result_folder = pathlib.Path(f\"results_largeGap\")\n",
+ "result_folder.mkdir(exist_ok=True)\n",
+ "for species_name, species in model_cur.sc.items:\n",
+ " results[species_name] = d.XDMFFile(\n",
+ " model_cur.mpi_comm_world, str(result_folder / f\"{species_name}.xdmf\")\n",
+ " )\n",
+ " results[species_name].parameters[\"flush_output\"] = True\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ "model_cur.to_pickle(\"model_cur.pkl\")\n",
+ "\n",
+ "# Set loglevel to warning in order not to pollute notebook output\n",
+ "logger.setLevel(logging.WARNING)\n",
+ "# Solve\n",
+ "displayed = False\n",
+ "while True:\n",
+ " # Solve the system\n",
+ " model_cur.monolithic_solve()\n",
+ " model_cur.adjust_dt()\n",
+ " # Save results for post processing\n",
+ " for species_name, species in model_cur.sc.items:\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ "\n",
+ " print(f\"Done with t={model_cur.t}\")\n",
+ " # End if we've passed the final time\n",
+ " if model_cur.t >= model_cur.final_t:\n",
+ " break\n",
+ "\n",
+ "# plt.plot(model_cur.tvec,)"
+ ]
+ }
+ ],
+ "metadata": {
+ "jupytext": {
+ "cell_metadata_filter": "-all",
+ "main_language": "python",
+ "notebook_metadata_filter": "-all"
+ },
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py
index e3b7335d..5eb54126 100644
--- a/smart/mesh_tools.py
+++ b/smart/mesh_tools.py
@@ -582,6 +582,215 @@ def meshSizeCallback(dim, tag, x, y, z, lc):
return (dmesh, mf2, mf3)
+def create_confined(
+ tubeRad: float = 10.0,
+ gapSize: float = 1.0,
+ cellVol: float = 1200.0,
+ hEdge: float = 0,
+ hInnerEdge: float = 0,
+ interface_marker: int = 12,
+ outer_marker: int = 10,
+ inner_vol_tag: int = 2,
+ outer_vol_tag: int = 1,
+ comm: MPI.Comm = d.MPI.comm_world,
+ verbose: bool = False,
+) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]:
+ """
+ Creates an axisymmetric mesh representing a cell within a cylindrical tube,
+ assuming axisymmetry about the r=0 axis.
+
+ Args:
+ tubeRad: Radius of tube
+ gapSize: Gap between cell and tube wall
+ cellVol: volume of cell (constrains overall geometry)
+ hEdge: maximum mesh size at the outer edge
+ hInnerEdge: maximum mesh size at the edge
+ of the inner compartment
+ interface_marker: The value to mark facets on the interface with
+ outer_marker: The value to mark facets on the outer ellipsoid with
+ inner_vol_tag: The value to mark the inner ellipsoidal volume with
+ outer_vol_tag: The value to mark the outer ellipsoidal volume with
+ comm: MPI communicator to create the mesh with
+ verbose: If true print gmsh output, else skip
+ Returns:
+ Tuple (mesh, facet_marker, cell_marker)
+ """
+ import gmsh
+
+ if tubeRad <= 0:
+ raise ValueError("Tube radius must be greater than 0")
+ if gapSize <= 0 or gapSize >= tubeRad:
+ raise ValueError("Gap between cell membrane and wall must be between 0 and tubeRad")
+ if cellVol <= 0:
+ raise ValueError("Cell volume must be greater than 0")
+
+ Lc = (cellVol - 4 * np.pi * (tubeRad - gapSize) ** 2) / (2 * np.pi * (tubeRad - gapSize))
+ zmax = Lc / 2 + (tubeRad - gapSize) + 10.0
+ rValsOuter = np.array([0.0, tubeRad, tubeRad, 0.0])
+ zValsOuter = np.array([zmax, zmax, -zmax, -zmax])
+ thetaVals = np.linspace(np.pi / 2, 0.0, 20)
+ rValsSphere = (tubeRad - gapSize) * np.cos(thetaVals)
+ zValsSphere = Lc / 2 + (tubeRad - gapSize) * np.sin(thetaVals)
+ maxOuterDim = zmax
+ maxInnerDim = Lc + (tubeRad - gapSize)
+
+ if np.isclose(hEdge, 0):
+ hEdge = 0.1 * maxOuterDim
+ if np.isclose(hInnerEdge, 0):
+ hInnerEdge = 0.2 * maxInnerDim
+ # Create the two axisymmetric body mesh using gmsh
+ gmsh.initialize()
+ gmsh.option.setNumber("General.Terminal", int(verbose))
+ gmsh.model.add("axisymm")
+ # first add outer body
+ outer_tag_list = []
+ outer_line_list = []
+ for i in range(len(rValsOuter)):
+ cur_tag = gmsh.model.occ.add_point(rValsOuter[i], 0, zValsOuter[i])
+ outer_tag_list.append(cur_tag)
+ if i > 0:
+ outer_line_list.append(gmsh.model.occ.add_line(cur_tag, outer_tag_list[-2]))
+ # include symm axis
+ outer_line_list.append(gmsh.model.occ.add_line(outer_tag_list[0], outer_tag_list[-1]))
+ outer_loop_tag = gmsh.model.occ.add_curve_loop(outer_line_list)
+ cell_plane_tag = gmsh.model.occ.add_plane_surface([outer_loop_tag])
+
+ # Add inner shape
+ inner_tag_list1 = []
+ for i in range(len(rValsSphere)):
+ cur_tag = gmsh.model.occ.add_point(rValsSphere[i], 0, zValsSphere[i])
+ inner_tag_list1.append(cur_tag)
+ inner_spline1_tag = gmsh.model.occ.add_spline(inner_tag_list1)
+ inner_tag_list2 = []
+ for i in range(len(rValsSphere)):
+ cur_tag = gmsh.model.occ.add_point(rValsSphere[-(i + 1)], 0, -zValsSphere[-(i + 1)])
+ inner_tag_list2.append(cur_tag)
+ inner_spline2_tag = gmsh.model.occ.add_spline(inner_tag_list2)
+ inner_cyl_line = gmsh.model.occ.add_line(inner_tag_list1[-1], inner_tag_list2[0])
+ symm_inner_tag = gmsh.model.occ.add_line(inner_tag_list1[0], inner_tag_list2[-1])
+ inner_loop_tag = gmsh.model.occ.add_curve_loop(
+ [inner_spline1_tag, inner_cyl_line, inner_spline2_tag, symm_inner_tag]
+ )
+ inner_plane_tag = gmsh.model.occ.add_plane_surface([inner_loop_tag])
+ cell_plane_list = [cell_plane_tag]
+ inner_plane_list = [inner_plane_tag]
+
+ outer_volume = []
+ inner_volume = []
+ all_volumes = []
+ inner_marker_list = []
+ outer_marker_list = []
+ for i in range(len(cell_plane_list)):
+ cell_plane_tag = cell_plane_list[i]
+ inner_plane_tag = inner_plane_list[i]
+ # Create interface between 2 objects
+ two_shapes, (outer_shape_map, inner_shape_map) = gmsh.model.occ.fragment(
+ [(2, cell_plane_tag)], [(2, inner_plane_tag)]
+ )
+ gmsh.model.occ.synchronize()
+
+ # Get the outer boundary
+ outer_shell = gmsh.model.getBoundary(two_shapes, oriented=False)
+ for i in range(len(outer_shell)):
+ outer_marker_list.append(outer_shell[i][1])
+ # Get the inner boundary
+ inner_shell = gmsh.model.getBoundary(inner_shape_map, oriented=False)
+ for i in range(len(inner_shell)):
+ inner_marker_list.append(inner_shell[i][1])
+ for tag in outer_shape_map:
+ all_volumes.append(tag[1])
+ for tag in inner_shape_map:
+ inner_volume.append(tag[1])
+
+ for vol in all_volumes:
+ if vol not in inner_volume:
+ outer_volume.append(vol)
+
+ # Add physical markers for facets
+ # set symmetry axis to 0 (no flux)
+ xmin, ymin, zmin = (-hInnerEdge / 10, -hInnerEdge / 10, -1)
+ xmax, ymax, zmax = (hInnerEdge / 10, hInnerEdge / 10, max(zValsOuter) + 1)
+ all_symm_bound = gmsh.model.occ.get_entities_in_bounding_box(
+ xmin, ymin, zmin, xmax, ymax, zmax, dim=1
+ )
+ symm_bound_markers = []
+ for i in range(len(all_symm_bound)):
+ symm_bound_markers.append(all_symm_bound[i][1])
+ # note that this first call sets the symmetry axis to tag 0 and
+ # this is not overwritten by the next calls to add_physical_group
+ gmsh.model.add_physical_group(1, symm_bound_markers, tag=0)
+ gmsh.model.add_physical_group(1, outer_marker_list, tag=outer_marker)
+ gmsh.model.add_physical_group(1, inner_marker_list, tag=interface_marker)
+
+ # Physical markers for "volumes"
+ gmsh.model.add_physical_group(2, outer_volume, tag=outer_vol_tag)
+ gmsh.model.add_physical_group(2, inner_volume, tag=inner_vol_tag)
+
+ def meshSizeCallback(dim, tag, x, y, z, lc):
+ # mesh length is hEdge at the PM and hInnerEdge at the inner membrane
+ # between these, the value is interpolated based on the relative distance
+ # between the two membranes.
+ # Inside the inner shape, the value is interpolated between hInnerEdge
+ # and lc3, where lc3 = max(hInnerEdge, 0.2*maxInnerDim)
+ # if innerRad=0, then the mesh length is interpolated between
+ # hEdge at the PM and 0.2*maxOuterDim in the center
+ lc1 = hEdge
+ lc2 = hInnerEdge
+ lc3 = max(hInnerEdge, 0.3 * (tubeRad - gapSize))
+ z_abs = np.abs(z)
+ if z_abs < Lc / 2:
+ if x > (tubeRad - gapSize):
+ in_outer = True
+ # dist_to_outer = tubeRad - x
+ dist_to_inner = x - (tubeRad - gapSize)
+ else:
+ in_outer = False
+ R_rel_inner = x / (tubeRad - gapSize)
+ else:
+ rTest = np.sqrt(x**2 + (z_abs - Lc / 2) ** 2)
+ if rTest > (tubeRad - gapSize):
+ in_outer = True
+ # dist_to_outer = tubeRad - x
+ dist_to_inner = rTest - (tubeRad - gapSize)
+ else:
+ in_outer = False
+ R_rel_inner = rTest / (tubeRad - gapSize)
+
+ if in_outer:
+ lcTest = lc2 + (lc1 - lc2) * (1 - np.exp(-dist_to_inner / 1.0))
+ else:
+ lcTest = lc2 + (lc3 - lc2) * (1 - R_rel_inner)
+ return lcTest
+
+ gmsh.model.mesh.setSizeCallback(meshSizeCallback)
+ # set off the other options for mesh size determination
+ gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
+ gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
+ gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
+ # this changes the algorithm from Frontal-Delaunay to Delaunay,
+ # which may provide better results when there are larger gradients in mesh size
+ gmsh.option.setNumber("Mesh.Algorithm", 5)
+
+ gmsh.model.mesh.generate(2)
+ rank = MPI.COMM_WORLD.rank
+ tmp_folder = pathlib.Path(f"tmp_2DCell_{rank}")
+ tmp_folder.mkdir(exist_ok=True)
+ gmsh_file = tmp_folder / "2DCell.msh"
+ gmsh.write(str(gmsh_file))
+ gmsh.finalize()
+
+ # return dolfin mesh of max dimension (parent mesh) and marker functions mf2 and mf3
+ dmesh, mf2, mf3 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 2, comm)
+ # ensure zero flux condition at r=0 axis
+ for f in d.facets(dmesh):
+ if np.isclose(f.midpoint().x(), 0.0):
+ mf2[f] = 0
+ # remove tmp mesh and tmp folder
+ gmsh_file.unlink(missing_ok=False)
+ tmp_folder.rmdir()
+ return (dmesh, mf2, mf3)
+
+
def create_cylinders(
outerRad: float = 1.0,
innerRad: float = 0.0,
diff --git a/smart/model.py b/smart/model.py
index 7964faff..2a25a1b1 100644
--- a/smart/model.py
+++ b/smart/model.py
@@ -255,6 +255,7 @@ def _init_1(self):
logger.debug("Checking validity of model (step 1 of ZZ)", extra=dict(format_type="title"))
self._init_1_1_check_mesh_dimensionality()
+ self._init_1_1b_CH_chem_potential_init()
self._init_1_2_check_namespace_conflicts()
self._init_1_3_check_parameter_dimensionality()
logger.debug(
@@ -372,6 +373,28 @@ def _init_1_1_check_mesh_dimensionality(self):
for compartment in self.cc:
compartment.is_volume = compartment.dimensionality == self.max_dim
+ def _init_1_1b_CH_chem_potential_init(self):
+ # all CH species require an additional species to be added to track chemical potential
+ # initial value is treated later
+ new_sp = []
+ for species in self.sc:
+ if species.CH:
+ init_chem_potential = 0.0
+ diff_units = self.cc[species.compartment_name].compartment_units ** 2 / unit.sec
+ chem_potential_sp = Species(
+ f"{species.name}_chem_potential",
+ init_chem_potential,
+ unit.dimensionless,
+ 0.0,
+ diff_units,
+ species.compartment_name,
+ )
+ chem_potential_sp.is_chem_potential = True
+ new_sp.append(chem_potential_sp)
+ species.chem_potential = chem_potential_sp
+ for sp in new_sp:
+ self.sc.add(sp)
+
def _init_1_2_check_namespace_conflicts(self):
"""Namespace checks:
@@ -584,6 +607,9 @@ def _init_2_4_check_for_unused_parameters_species_compartments(self):
all_parameters = set(chain.from_iterable([r.parameters for r in self.rc]))
all_species = set(chain.from_iterable([r.species for r in self.rc]))
+ for species in self.sc:
+ if species.is_chem_potential:
+ all_species.add(species.name)
all_compartments = set(chain.from_iterable([r.compartments for r in self.rc]))
if all_parameters != set(self.pc.keys):
print_str = (
@@ -987,6 +1013,9 @@ def _init_4_7_set_initial_conditions(self):
"""
logger.debug("Set function values to initial conditions", extra=dict(format_type="log"))
for species in self.sc:
+ if species.is_chem_potential:
+ species.D_dolfin = d.Constant(0.0) # set diffusion to zero (N/A)
+ continue # then initial condition is set to match concentration of assoc species
# first, initialize diffusion coefficient
if isinstance(species.D, float):
species.D_dolfin = d.Constant(species.D)
@@ -1024,6 +1053,21 @@ def _init_4_7_set_initial_conditions(self):
values[species.dof_map] = values_new[species.dof_map]
u_cur.vector().set_local(values)
u_cur.vector().apply("insert")
+ if species.CH:
+ lagrange = species.compartment.deform_logic or species.alt_deform_logic
+ if lagrange:
+ logger.error("CH species are not compatible with Lagrange approach yet!")
+ for ckey in species.chem_potential.u.keys():
+ A_hat = species.A_hat
+ phi_cur = species.u[ckey] / species.umax
+ cfunc = (
+ d.ln(phi_cur)
+ - d.ln(1 - phi_cur)
+ - A_hat * (2 * phi_cur - 1)
+ + (A_hat / species.umax) * d.div(d.grad(phi_cur))
+ )
+ Vc = species.chem_potential.V
+ species.chem_potential.u[ckey].assign(d.project(cfunc, Vc))
species.alt_vel_logic = np.any([vel != 0.0 for vel in species.alt_vel])
species.alt_deform_logic = np.any([deform != 0.0 for deform in species.alt_deform])
if species.alt_vel_logic and species.alt_deform_logic:
@@ -1339,8 +1383,53 @@ def _init_5_2_create_variational_forms(self):
J = d.Constant(1.0)
if self.config.flags["axisymmetric_model"]:
J = x[0] * J
+ # catch CH case
+ if species.CH:
+ if species.is_chem_potential:
+ logger.debug(
+ "Chemical potential equation is defined with concentration,"
+ "skipping to next species"
+ )
+ continue
+ else:
+ u_c = species.chem_potential._usplit["u"]
+ v_c = species.chem_potential.v
+ A_hat = species.A_hat
+ phi_cur = u / species.umax
+ df_c = d.ln(phi_cur) - d.ln(1 - phi_cur) - A_hat * (2 * phi_cur - 1)
+ # CForm scaling factor
+ CScale = (float(species.D) * species.umax) / (4 * np.pi) # assuming l^2 = 4*pi
+ CForm = J * (
+ (u_c - df_c) * v_c * dx
+ - (A_hat / species.umax) * d.inner(d.grad(phi_cur), d.grad(v_c)) * dx
+ )
+ # chemical potential is in units of kBT for convenience
+ Dform = J * D * u * d.inner(d.grad(u_c), d.grad(v)) * dx
+ self.forms.add(
+ Form(
+ f"chem_potential_{species.name}",
+ CForm,
+ species.chem_potential,
+ "chem_potential",
+ Dform_units,
+ True,
+ linear_wrt_comp,
+ form_scaling=CScale,
+ )
+ )
+ self.forms.add(
+ Form(
+ f"diffusion_{species.name}",
+ Dform,
+ species,
+ "diffusion",
+ Dform_units,
+ True,
+ linear_wrt_comp,
+ )
+ )
# diffusion term
- if species.D == 0:
+ elif species.D == 0:
logger.debug(
f"Species {species.name} has a diffusion coefficient of 0. "
"Skipping creation of diffusive form.",
@@ -1368,8 +1457,6 @@ def _init_5_2_create_variational_forms(self):
)
else:
Dform = J * D * d.inner(d.grad(u), d.grad(v)) * dx
- # exponent is -2 because of two gradients
-
self.forms.add(
Form(
f"diffusion_{species.name}",
diff --git a/smart/model_assembly.py b/smart/model_assembly.py
index 55abd2b2..ffa79e54 100644
--- a/smart/model_assembly.py
+++ b/smart/model_assembly.py
@@ -988,6 +988,11 @@ class Species(ObjectInstance):
alt_deform: list = dataclasses.field(default_factory=lambda: [0.0, 0.0, 0.0])
alt_vel: list = dataclasses.field(default_factory=lambda: [0.0, 0.0, 0.0])
alt_manual_update: bool = False
+ CH: bool = False
+ is_chem_potential: bool = False
+ A_hat: float = 0.0
+ umax: float = 0.0
+ # if this is a CH conc, then we also need fields: chem_potential, A_hat, umax
def to_dict(self):
"Convert to a dict that can be used to recreate the object."
@@ -1061,6 +1066,12 @@ def __post_init__(self):
else:
raise TypeError("Diffusion coefficient must a float, int, or string")
+ if self.CH:
+ if not hasattr(self, "A_hat"):
+ raise ValueError("A_hat must be provided for CH species")
+ if not hasattr(self, "umax"):
+ raise ValueError("umax must be provided for CH variable")
+
self._convert_pint_quantity_to_unit()
self._check_input_type_validity()
self._convert_pint_unit_to_quantity()