diff --git a/examples/example2/example2.ipynb b/examples/example2/example2.ipynb index 46b7a3c..5555dfa 100644 --- a/examples/example2/example2.ipynb +++ b/examples/example2/example2.ipynb @@ -429,6 +429,7 @@ " outer_tag=surf_tag, outer_marker=edge_tag)\n", " mesh_folder = pathlib.Path(f\"ellipse_mesh_AR{aspect_ratios[i]}\")\n", " mesh_folder.mkdir(exist_ok=True)\n", + " mesh_file = mesh_folder / \"ellipse_mesh.h5\"\n", " mesh_tools.write_mesh(ellipse_mesh, mf1, mf2, mesh_file)\n", " parent_mesh = mesh.ParentMesh(\n", " mesh_filename=str(mesh_file),\n", diff --git a/examples/example2/example2_lagrangeTest.py b/examples/example2/example2_lagrangeTest.py new file mode 100644 index 0000000..28e05d7 --- /dev/null +++ b/examples/example2/example2_lagrangeTest.py @@ -0,0 +1,257 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Example 2: Simple 2D cell signaling model +# Test alterations in mesh vs. Lagrangian mapping for different cases + +import dolfin as d +import numpy as np +import pathlib +import logging + +from smart import config, mesh, model, mesh_tools +from smart.units import unit +from smart.model_assembly import ( + Compartment, + Parameter, + Reaction, + Species, + SpeciesContainer, + ParameterContainer, + CompartmentContainer, + ReactionContainer, +) +from matplotlib import pyplot as plt + +logger = logging.getLogger("smart") +logger.setLevel(logging.INFO) + +um = unit.um +molecule = unit.molecule +sec = unit.sec +dimensionless = unit.dimensionless +D_unit = um**2 / sec +surf_unit = molecule / um**2 +flux_unit = molecule / (um * sec) +edge_unit = molecule / um + +# ============================================================================================= +# Compartments +# ============================================================================================= +# name, topological dimensionality, length scale units, marker value +Cyto = Compartment("Cyto", 2, um, 1) +PM = Compartment("PM", 1, um, 3) +cc = CompartmentContainer() +cc.add([Cyto, PM]) + +# ============================================================================================= +# Species +# ============================================================================================= +# name, initial concentration, concentration units, diffusion, diffusion units, compartment +A = Species("A", 1.0, surf_unit, 0.001, D_unit, "Cyto") +X = Species("X", 1.0, edge_unit, 0.001, D_unit, "PM") +B = Species("B", 0.0, edge_unit, 0.001, D_unit, "PM") +sc = SpeciesContainer() +sc.add([A, X, B]) + +# ============================================================================================= +# Parameters and Reactions +# ============================================================================================= + +# Reaction of A and X to make B (Cyto-PM reaction) +kon = Parameter("kon", 1.0, 1 / (surf_unit * sec)) +koff = Parameter("koff", 1.0, 1 / sec) +r1 = Reaction( + "r1", + ["A", "X"], + ["B"], + param_map={"on": "kon", "off": "koff"}, + species_map={"A": "A", "X": "X", "B": "B"}, +) + +pc = ParameterContainer() +pc.add([kon, koff]) +rc = ReactionContainer() +rc.add([r1]) + +h_ellipse = 0.2 +xrad = 1.0 +yrad = 1.0 +surf_tag = 1 +edge_tag = 3 +config_cur = config.Config() +config_cur.solver.update( + { + "final_t": 2.0, + "initial_dt": 0.1, + "time_precision": 6, + } +) + +# iterate over additional aspect ratios, enforcing the same ellipsoid area in all cases +aspect_ratios = [1, 1.5**2, 4, 9, 16] +l2_lagrange = [] +l2_eul = [] +Avecs = [] + +parent_folder = pathlib.Path("lagrange_testing") +parent_folder.mkdir(exist_ok=True) + +for i in range(len(aspect_ratios)): + + udefs = [] + lagrangeSS = [] + refSS = [] + for lagrange in [False, True]: + if lagrange: + Cyto = Compartment( + "Cyto", + 2, + um, + 1, + deform=( + f"x*{np.sqrt(aspect_ratios[i])-1}", + f"y*{1/np.sqrt(aspect_ratios[i]) - 1}", + "0", + ), + ) + PM = Compartment( + "PM", + 1, + um, + 3, + deform=( + f"x*{np.sqrt(aspect_ratios[i])-1}", + f"y*{1/np.sqrt(aspect_ratios[i]) - 1}", + "0", + ), + ) + else: + Cyto = Compartment("Cyto", 2, um, 1) + PM = Compartment("PM", 1, um, 3) + cc = CompartmentContainer() + cc.add([Cyto, PM]) + + # Create mesh + if not lagrange: + xrad = 1.0 * np.sqrt(aspect_ratios[i]) + yrad = 1.0 / np.sqrt(aspect_ratios[i]) + ellipse_mesh, mf1, mf2 = mesh_tools.create_ellipses( + xrad, + yrad, + hEdge=h_ellipse, + hInnerEdge=h_ellipse, + outer_tag=surf_tag, + outer_marker=edge_tag, + ) + mesh_folder = parent_folder / f"ellipse_mesh_AR{aspect_ratios[i]}" + mesh_folder.mkdir(exist_ok=True) + else: + mesh_folder = parent_folder / "ellipse_mesh_AR1" + + mesh_refinements = [0, 1, 2, 3, 4, 5, 6] + for j in range(len(mesh_refinements)): + mesh_file = mesh_folder / f"ellipse_mesh_{mesh_refinements[j]}.h5" + if not lagrange: + if mesh_refinements[j] > 0: + d.parameters["refinement_algorithm"] = "plaza_with_parent_facets" + ellipse_mesh = d.adapt(ellipse_mesh) + mf2 = d.adapt(mf2, ellipse_mesh) + mf1 = d.adapt(mf1, ellipse_mesh) + mesh_tools.write_mesh(ellipse_mesh, mf1, mf2, mesh_file) + + parent_mesh = mesh.ParentMesh( + mesh_filename=str(mesh_file), + mesh_filetype="hdf5", + name="parent_mesh", + ) + # init model + model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh) + model_cur.initialize() + results = dict() + result_folder = ( + parent_folder + / f"resultsEllipse_AR{aspect_ratios[i]}_lagrange{lagrange}_{mesh_refinements[j]}" + ) + result_folder.mkdir(exist_ok=True) + for species_name, species in model_cur.sc.items: + results[species_name] = d.XDMFFile( + model_cur.mpi_comm_world, str(result_folder / f"{species_name}.xdmf") + ) + results[species_name].parameters["flush_output"] = True + results[species_name].write(model_cur.sc[species_name].u["u"], model_cur.t) + if lagrange: + u1file = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / "u1var.xdmf")) + u1file.parameters["flush_output"] = True + u1file.write(cc["Cyto"].deform_func, model_cur.t) + udef = cc["Cyto"].deform_func + Fcur = d.Identity(3) + d.grad(udef) + Jcur = d.det(Fcur) + else: + Jcur = d.Constant(1.0) + avg_A = [A.initial_condition] + dx = d.Measure("dx", domain=model_cur.cc["Cyto"].dolfin_mesh) + volume = d.assemble_mixed(Jcur * dx) + while True: + # Solve the system + model_cur.monolithic_solve() + # Save results for post processing + for species_name, species in model_cur.sc.items: + results[species_name].write(model_cur.sc[species_name].u["u"], model_cur.t) + if lagrange: + u1file.write(cc["Cyto"].deform_func, model_cur.t) + int_val = d.assemble_mixed(Jcur * model_cur.sc["A"].u["u"] * dx) + avg_A.append(int_val / volume) + # End if we've passed the final time + if model_cur.t >= model_cur.final_t: + break + + np.savetxt(str(result_folder / "avg_A.txt"), avg_A) + Avecs.append(avg_A) + + if lagrange: + plt.plot( + model_cur.tvec, + avg_A, + linestyle="dashed", + label=f"Aspect ratio = {aspect_ratios[i]}, Lagrange", + ) + lagrangeSS.append(model_cur.sc["A"].sol.copy()) + udefs.append(udef) + else: + plt.plot(model_cur.tvec, avg_A, label=f"Aspect ratio = {aspect_ratios[i]}") + refSS.append(model_cur.sc["A"].sol.copy()) + + refFcn = refSS[-1] # use the most refined reference (eulerian) mesh + Vcur = refFcn.function_space() + coords = Vcur.tabulate_dof_coordinates() + dx = d.Measure("dx", domain=Vcur.mesh()) + l2_lagrangeCur = [] + l2_refCur = [] + for k in range(len(refSS) - 1): + lagrange_fun = d.Function(Vcur) + lagrange_vec = lagrange_fun.vector().get_local() + eul_fun = d.Function(Vcur) + eul_vec = eul_fun.vector().get_local() + + refSS[k].set_allow_extrapolation(True) + Lmesh = lagrangeSS[k].function_space().mesh() + d.ALE.move(Lmesh, udefs[k]) + Lmesh.bounding_box_tree().build(Lmesh) + lagrangeSS[k].set_allow_extrapolation(True) + for idx in range(len(lagrange_vec)): + cur_coord = coords[idx] + eul_vec[idx] = refSS[k](cur_coord) + lagrange_vec[idx] = lagrangeSS[k](cur_coord) + eul_fun.vector().set_local(eul_vec) + eul_fun.vector().apply("insert") + lagrange_fun.vector().set_local(lagrange_vec) + lagrange_fun.vector().apply("insert") + l2_lagrangeCur.append(np.sqrt(d.assemble((lagrange_fun - refFcn) ** 2 * dx))) + l2_refCur.append(np.sqrt(d.assemble((eul_fun - refFcn) ** 2 * dx))) + + np.savetxt(str(parent_folder / f"l2_lagrange_AR{aspect_ratios[i]}.txt"), l2_lagrangeCur) + np.savetxt(str(parent_folder / f"l2_eul_AR{aspect_ratios[i]}.txt"), l2_refCur) + +plt.legend() +plt.show() diff --git a/examples/example2/example2_withadvection.ipynb b/examples/example2/example2_withadvection.ipynb new file mode 100644 index 0000000..1944777 --- /dev/null +++ b/examples/example2/example2_withadvection.ipynb @@ -0,0 +1,498 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "f65f18d7", + "metadata": {}, + "source": [ + "# Example 2: Simple 2D cell signaling model\n", + "\n", + "We model a reaction between the cell interior and cell membrane in a 2D geometry:\n", + "- Cyto - 2D cell \"volume\"\n", + "- PM - 1D cell boundary (represents plasma membrane)\n", + "\n", + "Model from [Rangamani et al, 2013, Cell](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3874130/). A cytosolic species, \"A\", reacts with a species on the PM, \"B\", to form a new species on the PM, \"X\". The resulting PDE and boundary condition for species A are as follows:\n", + "\n", + "$$\n", + "\\frac{\\partial{C_A}}{\\partial{t}} = D_A \\nabla ^2 C_A \\quad \\text{in} \\; \\Omega_{Cyto}\\\\\n", + "\\text{B.C. for A:} \\quad D_A (\\textbf{n} \\cdot \\nabla C_A) = -k_{on} C_A N_X + k_{off} N_B \\quad \\text{on} \\; \\Gamma_{PM}\n", + "$$\n", + "\n", + "Similarly, the PDEs for X and B are given by:\n", + "$$\n", + "\\frac{\\partial{N_X}}{\\partial{t}} = D_X \\nabla ^2 N_X - k_{on} C_A N_X + k_{off} N_B \\quad \\text{on} \\; \\Gamma_{PM}\\\\\n", + "\\frac{\\partial{N_B}}{\\partial{t}} = D_B \\nabla ^2 N_B + k_{on} C_A N_X - k_{off} N_B \\quad \\text{on} \\; \\Gamma_{PM}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b224bea7", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "import matplotlib.image as mpimg\n", + "img_A = mpimg.imread('axb-diagram.png')\n", + "plt.imshow(img_A)\n", + "plt.axis('off')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "a59f3428", + "metadata": {}, + "source": [ + "Imports and logger initialization:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc398816", + "metadata": { + "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, 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", + " Form,\n", + ")\n", + "from matplotlib import pyplot as plt\n", + "import matplotlib.image as mpimg\n", + "\n", + "logger = logging.getLogger(\"smart\")\n", + "logger.setLevel(logging.INFO)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "95b9d865", + "metadata": {}, + "source": [ + "First, we define the various units for use in the model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f4023cf", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "um = unit.um\n", + "molecule = unit.molecule\n", + "sec = unit.sec\n", + "dimensionless = unit.dimensionless\n", + "D_unit = um**2 / sec\n", + "surf_unit = molecule / um**2\n", + "flux_unit = molecule / (um * sec)\n", + "edge_unit = molecule / um" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "46582d26", + "metadata": {}, + "source": [ + "Next we generate the model by assembling the compartment, species, parameter, and reaction containers (see Example 1 or API documentation for more details)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "09079b17", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "# =============================================================================================\n", + "# Compartments\n", + "# =============================================================================================\n", + "# name, topological dimensionality, length scale units, marker value\n", + "# Cyto = Compartment(\"Cyto\", 2, um, 1, vel=(\"-0.1*y*(1 - (x**2+y**2))\",\"0.1*x*(1 - (x**2+y**2))\",\"0\"))\n", + "# Cyto = Compartment(\"Cyto\", 2, um, 1, deform=(\"-0.1*y*t*(1 - (x**2+y**2))\",\"0.1*x*t*(1 - (x**2+y**2))\",\"0\"), manual_update=True)\n", + "Cyto = Compartment(\"Cyto\", 2, um, 1, deform=(\"((1+sign(t-10))/2)*0.005*(t-10)*x/(1-0.005*(t-10))\",\n", + " \"-((1+sign(t-10))/2)*0.005*(t-10)*y\",\"0\")) # dilation\n", + "PM = Compartment(\"PM\", 1, um, 3, deform=(\"((1+sign(t-10))/2)*0.005*(t-10)*x/(1-0.005*(t-10))\",\n", + " \"-((1+sign(t-10))/2)*0.005*(t-10)*y\",\"0\"))\n", + "cc = CompartmentContainer()\n", + "cc.add([Cyto, PM])\n", + "\n", + "# =============================================================================================\n", + "# Species\n", + "# =============================================================================================\n", + "# name, initial concentration, concentration units, diffusion, diffusion units, compartment\n", + "# A = Species(\"A\", \"10*exp(-(x**2+(y-0.8)**2)/0.2**2)\", surf_unit, 0.001, D_unit, \"Cyto\")\n", + "A = Species(\"A\", 1.0, surf_unit, 0.001, D_unit, \"Cyto\")\n", + "X = Species(\"X\", 1.0, edge_unit, 0.001, D_unit, \"PM\")\n", + "B = Species(\"B\", 0.0, edge_unit, 0.001, D_unit, \"PM\")\n", + "sc = SpeciesContainer()\n", + "sc.add([A, X, B])\n", + "\n", + "# =============================================================================================\n", + "# Parameters and Reactions\n", + "# =============================================================================================\n", + "\n", + "# Reaction of A and X to make B (Cyto-PM reaction)\n", + "kon = Parameter(\"kon\", 1.0, 1/(surf_unit*sec))\n", + "koff = Parameter(\"koff\", 1.0, 1/sec)\n", + "r1 = Reaction(\"r1\", [\"A\", \"X\"], [\"B\"],\n", + " param_map={\"on\": \"kon\", \"off\": \"koff\"},\n", + " species_map={\"A\": \"A\", \"X\": \"X\", \"B\": \"B\"})\n", + "\n", + "pc = ParameterContainer()\n", + "pc.add([kon, koff])\n", + "rc = ReactionContainer()\n", + "rc.add([r1])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "15c35d39", + "metadata": {}, + "source": [ + "Now we create a circular mesh (mesh built using gmsh in `smart.mesh_tools`), along with marker functions `mf2` and `mf1`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe56e162", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "# Create mesh\n", + "h_ellipse = 0.1\n", + "xrad = 1.0\n", + "yrad = 1.0\n", + "surf_tag = 1\n", + "edge_tag = 3\n", + "ellipse_mesh, mf1, mf2 = mesh_tools.create_ellipses(2*xrad, 0.5*yrad, hEdge=h_ellipse,\n", + " outer_tag=surf_tag, outer_marker=edge_tag)\n", + "visualization.plot_dolfin_mesh(ellipse_mesh, mf2, view_xy=True)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "fe04ad6b", + "metadata": {}, + "source": [ + "Write mesh and meshfunctions to file, then create `mesh.ParentMesh` object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e15255a1", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "mesh_folder = pathlib.Path(\"ellipse_mesh_AR1\")\n", + "mesh_folder.mkdir(exist_ok=True)\n", + "mesh_file = mesh_folder / \"ellipse_mesh.h5\"\n", + "mesh_tools.write_mesh(ellipse_mesh, mf1, mf2, mesh_file)\n", + "\n", + "parent_mesh = mesh.ParentMesh(\n", + " mesh_filename=str(mesh_file),\n", + " mesh_filetype=\"hdf5\",\n", + " name=\"parent_mesh\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d1c0cab2", + "metadata": {}, + "source": [ + "Initialize model and solvers." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b059df37", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "config_cur = config.Config()\n", + "config_cur.solver.update(\n", + " {\n", + " \"final_t\": 130.0,\n", + " \"initial_dt\": 1.0,\n", + " \"time_precision\": 6,\n", + " }\n", + ")\n", + "\n", + "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n", + "model_cur.initialize()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "e97c54a3", + "metadata": {}, + "source": [ + "Save model information to .pkl file and write initial conditions to file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "274cc41d", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "model_cur.to_pickle('model_cur.pkl')\n", + "results = dict()\n", + "result_folder = pathlib.Path(\"resultsEllipse_testShapeChange\")\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", + "# save deformation or velocity field as well\n", + "if cc[\"Cyto\"].vel_logic:\n", + " velfile = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"vel.xdmf\"))\n", + " velfile.parameters[\"flush_output\"] = True\n", + " velfile.write(sc[\"A\"].compartment.vel_func, model_cur.t)\n", + "elif cc[\"Cyto\"].deform_logic:\n", + " u1file = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"u1var.xdmf\"))\n", + " u1file.parameters[\"flush_output\"] = True\n", + " u1file.write(cc[\"Cyto\"].deform_func, model_cur.t)\n", + " u2file = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"u2var.xdmf\"))\n", + " u2file.parameters[\"flush_output\"] = True\n", + " u2file.write(cc[\"PM\"].deform_func, model_cur.t)\n", + " # v1file = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"v1.xdmf\"))\n", + " # v1file.parameters[\"flush_output\"] = True\n", + " # v2file = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"v2.xdmf\"))\n", + " # v2file.parameters[\"flush_output\"] = True" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "355a426e", + "metadata": {}, + "source": [ + "Solve the system until `model_cur.t > model_cur.final_t`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40b213ea", + "metadata": { + "metadata": {} + }, + "outputs": [], + "source": [ + "tvec = [0]\n", + "dx_Cyto = d.Measure(\"dx\", domain=model_cur.cc['Cyto'].dolfin_mesh)\n", + "dx_PM = d.Measure(\"dx\", domain=model_cur.cc['PM'].dolfin_mesh)\n", + "volume = d.assemble_mixed(1.0*dx_Cyto)\n", + "PM_SA = d.assemble_mixed(1.0*dx_PM)\n", + "avg_A = [d.assemble_mixed(A.u['u']*dx_Cyto) / volume]\n", + "avg_X = [d.assemble_mixed(X.u[\"u\"]*dx_PM) / PM_SA]\n", + "avg_B = [d.assemble_mixed(B.u[\"u\"]*dx_PM) / PM_SA]\n", + "# Set loglevel to warning in order not to pollute notebook output\n", + "logger.setLevel(logging.WARNING)\n", + "x = d.SpatialCoordinate(sc[\"A\"].compartment.dolfin_mesh)\n", + "# if cc[\"Cyto\"].deform_logic:\n", + "# udef = cc[\"Cyto\"].deform_func\n", + "# udef_local = udef.copy()\n", + "# V_vector = udef.function_space()\n", + "# vel_expr = d.Expression((\"-0.1*x[1]*(1 - (pow(x[0],2)+pow(x[1],2)))\",\n", + "# \"0.1*x[0]*(1 - (pow(x[0],2)+pow(x[1],2)))\",\"0\"), degree=1)\n", + "# vel = d.interpolate(vel_expr, V_vector)\n", + " # v1file.write(vel, model_cur.t)\n", + " # vel_approx = d.project((udef - cc[\"Cyto\"].deform_prev) / model_cur.dT, V_vector)\n", + " # v2file.write(vel_approx, model_cur.t)\n", + "\n", + "while True:\n", + " # Solve the system\n", + " model_cur.monolithic_solve()\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", + " avg_A.append(d.assemble_mixed(A.u['u']*dx_Cyto) / volume)\n", + " avg_X.append(d.assemble_mixed(X.u[\"u\"]*dx_PM) / PM_SA)\n", + " avg_B.append(d.assemble_mixed(B.u[\"u\"]*dx_PM) / PM_SA)\n", + " tvec.append(model_cur.t)\n", + " # save deformation or velocity\n", + " if cc[\"Cyto\"].vel_logic:\n", + " velfile.write(sc[\"A\"].compartment.vel_func, model_cur.t)\n", + " elif cc[\"Cyto\"].deform_logic:\n", + " # cc[\"Cyto\"].deform_prev.assign(udef_local)\n", + " # udef.assign(udef_local + vel*d.Constant(model_cur.tvec[-1]-model_cur.tvec[-2]))\n", + " u1file.write(cc[\"Cyto\"].deform_func, model_cur.t)\n", + " u2file.write(cc[\"PM\"].deform_func, model_cur.t)\n", + " # multCur = d.project(model_cur.fc[\"r1 [X (f)]\"].integral_factor,X.V)\n", + " # if np.any(np.isnan(multCur.vector()[:])):\n", + " # udef = cc[\"Cyto\"].deform_func\n", + " # Fcur = d.Identity(3) + d.grad(udef)\n", + " # Jcur = d.det(Fcur)\n", + " # Nexpr = d.Expression((\"x[0]/R\", \"x[1]/R\", \"0.0\"), degree=1, R=1)\n", + " # test_factor = Jcur*d.sqrt(d.inner(d.dot(Nexpr,d.inv(Fcur)),d.dot(Nexpr,d.inv(Fcur))))\n", + " # Vcur = d.VectorFunctionSpace(X.compartment.dolfin_mesh, \"P\", 1)\n", + " # N = d.interpolate(Nexpr, Vcur)\n", + " # N = self.surface.normals\n", + " \n", + " # for name, flux in model_cur.fc.items:\n", + " # if \"r1\" in name:\n", + " # Vcur = flux.integral_factor.function_space()\n", + " # not_assigned = True\n", + " # while not_assigned:\n", + " # try:\n", + " # flux.integral_factor.assign(d.project(test_factor, Vcur))\n", + " # not_assigned = False\n", + " # except:\n", + " # print('Try again!')\n", + " # multfile.write(model_cur.fc[\"r1 [X (f)]\"].integral_factor, model_cur.t)\n", + " # udef_local = udef.copy()\n", + " # # compute new velocities\n", + " # xcur = x[0] + udef[0]\n", + " # ycur = x[1] + udef[1]\n", + " # vcur_x = d.project(-0.1*ycur*(1 - (xcur**2+ycur**2)), V_vector.sub(0).collapse())\n", + " # vcur_y = d.project(0.1*xcur*(1 - (xcur**2+ycur**2)), V_vector.sub(1).collapse())\n", + " # vel.vector()[V_vector.sub(0).dofmap().dofs()] = vcur_x.vector()\n", + " # vel.vector()[V_vector.sub(1).dofmap().dofs()] = vcur_y.vector()\n", + " # vel.vector().apply(\"insert\")\n", + " # v1file.write(vel, model_cur.t)\n", + " # vel_approx.assign(d.project((udef - cc[\"Cyto\"].deform_prev) / model_cur.dT, V_vector))\n", + " # v2file.write(vel_approx, 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" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "52c0fd7e", + "metadata": {}, + "source": [ + "Now we plot the concentration of A in the cell over time and compare this to analytical predictions for a high value of the diffusion coefficient. As $D_A \\rightarrow \\infty$, the steady state concentration of A will be given by the positive root of the following polynomial:\n", + "\n", + "$$\n", + "-k_{on} c_A ^2 - \\left( k_{on} c_{X} (t=0) \\frac{SA_{PM}}{vol_{cyto}} + k_{off} - k_{on} c_A (t=0) \\right) c_A + k_{off} c_A(t=0)\n", + "$$\n", + "\n", + "Note that in this 2D case, $SA_{PM}$ is the perimeter of the ellipse and $vol_{cyto}$ is the area of the ellipse. These can be thought of as a surface area and volume if we extrude the 2D shape by some characteristic thickness." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4afd3149", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(tvec, avg_A, label='SMART simulation')\n", + "plt.xlabel('Time (s)')\n", + "plt.ylabel('A concentration $\\mathrm{(molecules/μm^2)}$')\n", + "SA_vol = 4/(xrad + yrad)\n", + "root_vals = np.roots([-kon.value,\n", + " -kon.value*avg_X[0]*SA_vol - koff.value + kon.value*avg_A[0],\n", + " koff.value*avg_A[0]])\n", + "ss_pred = root_vals[root_vals > 0]\n", + "plt.plot(tvec, np.ones(len(avg_A))*ss_pred, '--', label='Steady-state analytical prediction')\n", + "plt.legend()\n", + "# from cowpy import cow\n", + "# if True:#percent_error_analytical > 0.1:\n", + "# stegy = cow.Stegosaurus(thoughts=True)\n", + "# stegy_msg = stegy.milk(\"Failed test: Example 2 results deviate from the analytical prediction\")\n", + "# print(stegy_msg)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "05f30a3c", + "metadata": {}, + "source": [ + "Plot A concentration in the cell at the final time point." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffae481c", + "metadata": {}, + "outputs": [], + "source": [ + "visualization.plot(model_cur.sc[\"A\"].u[\"u\"], view_xy=True)" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/example5/example5_withSOCE.ipynb b/examples/example5/example5_withSOCE.ipynb new file mode 100644 index 0000000..198c844 --- /dev/null +++ b/examples/example5/example5_withSOCE.ipynb @@ -0,0 +1,854 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "f65f18d7", + "metadata": {}, + "source": [ + "# Example 5: Generic cell signaling system in 3D, illustrating a SOCE-like effect\n", + "\n", + "This example is slightly altered from the normal version of example 5 to illustrate the role of store-operated calcium entry (SOCE) accompanying calcium release from an interior compartment.\n", + "\n", + "Geometry is divided into 4 domains; two volumes, and two surfaces:\n", + "- cytosol (Cyto): $\\Omega_{Cyto}$\n", + "- endoplasmic reticulum volume (ER): $\\Omega_{ER}$\n", + "- plasma membrane (PM): $\\Gamma_{PM}$\n", + "- ER membrane (ERm): $\\Gamma_{ERm}$\n", + "\n", + "For simplicity, here we consider a \"cube-within-a-cube\" geometry, in which the smaller\n", + "inner cube represents a section of ER and one face of the outer cube ($x=0$) represents the PM. The other\n", + "faces of the outer cube are treated as no flux boundaries. The space outside\n", + "the inner cube but inside the outer cube is classified as cytosol.\n", + "\n", + "There are three function-spaces on these three domains:\n", + "\n", + "$$\n", + "u^{Cyto} = [A, B] \\quad \\text{on} \\quad \\Omega^{Cyto}\\\\\n", + "u^{ER} = [AER] \\quad \\text{on} \\quad \\Omega^{ER}\\\\\n", + "v^{ERm} = [R, Ro] \\quad \\text{on} \\quad \\Gamma^{ERm}\n", + "$$\n", + "\n", + "In words, this says that species A and B reside in the cytosolic volume, \n", + "species AER corresponds to an amount of species A that lives in the ER volume,\n", + "and species R (closed receptor/channel) and Ro (open receptor/channel) reside on the ER membrane.\n", + "\n", + "In this model, species B reacts with a receptor/channel, R, on the ER membrane, causing it to open (change state from R->Ro), \n", + "allowing species A to flow out of the ER and into the cytosol. \n", + "Note that this is roughly similar to an IP3 pulse at the PM, leading to Ca2+ release from the ER,\n", + "where, by analogy, species B is similar to IP3 and species A is similar to Ca2+. A more comprehensive\n", + "model of Ca2+ dynamics in particular is implemented in Example 6.\n", + "\n", + "Here, we also introduce a flux analogous to SOCE, wherein the flux of A through the PM depends on depletion of A from the ER.\n", + "For simplicity, we scale this flux with the distance between the PM and ER ($d_{PM-ER}$) and the concentration of A at closest point of ER to the PM ($A_{ER,close}$):\n", + "\n", + "$$\n", + "J_{SOCE} = \\frac{d_0}{d_{PM-ER} + d_0} \\frac{1}{1 + \\left(\\frac{A_{ER}}{A_{ref}}\\right)^4},\n", + "$$\n", + "\n", + "where $d_0$ is a characteristic length scale controlling the likelihood of STIM-Orai binding and $A_{ref}$ governs the ER calcium depletion required to trigger SOCE." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5edff5d", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "import matplotlib.image as mpimg\n", + "img_A = mpimg.imread('example5-diagram.png')\n", + "plt.imshow(img_A)\n", + "plt.axis('off')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "f543d75e", + "metadata": {}, + "source": [ + "As specified in our [mathematical documentation](https://rangamanilabucsd.github.io/smart/docs/math.html), assuming diffusive transport, the PDE and boundary condition for each of these volumetric species takes the form:\n", + "\n", + "$$\n", + "\\partial_t u_i^m - \\nabla \\cdot ( D_i^m \\nabla (u_i^m) ) - f_i^m(u_i^m) = 0 \\qquad \\text{ in } \\Omega^m\\\\\n", + "D_i \\nabla u_i^m \\cdot n^m - R_i^{q} (u^m, u^n, v^q) = 0 \\qquad \\text{ on } \\Gamma^{q}\n", + "$$\n", + "\n", + "and the surface species take the form:\n", + "\n", + "$$\n", + "\\partial_t v_i^q - \\nabla_S \\cdot (D_i^q \\nabla_S v ) - g_i^q ( u^m, u^n, v^q ) = 0 \\qquad \\text{ on } \\Gamma^{q}\\\\\n", + "D_i \\nabla v_i^q \\cdot n^q = 0 \\qquad \\text{ on } \\partial\\Gamma^{q}\n", + "$$\n", + "\n", + "Our reaction terms and boundary conditions are chosen according to the system described above. For the purposes of this simplified example we use linear mass action in all reaction terms except SOCE. Explicitly writing out this system of PDEs, we have:\n", + "\n", + "\\begin{align}\n", + " \\partial_t u_B^{Cyto} - D_B^{Cyto} \\nabla^2 u_B^{Cyto} + k_{2f} u_B^{Cyto} &= 0 \\qquad \\text{ in } \\Omega^{Cyto}\\\\\n", + " D_B^{Cyto} \\nabla u_B^{Cyto} \\cdot n^{Cyto} + j_1[t] &= 0 \\qquad \\text{ on } \\Gamma^{PM} \\nonumber\\\\\n", + " D_B^{Cyto} \\nabla u_B^{Cyto} \\cdot n^{Cyto} + k_{3f} v_R^{ERm} u_B^{Cyto} - k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t u_A^{Cyto} - D_A^{Cyto} \\nabla^2 u_A^{Cyto} &= 0 \\qquad \\text{ in } \\Omega^{Cyto}\\\\\n", + " D_A^{Cyto} \\nabla u_A^{Cyto} \\cdot n^{Cyto} - J_{SOCE} &= 0 \\qquad \\text{ on } \\Gamma^{PM} \\nonumber\\\\\n", + " D_A^{Cyto} \\nabla u_A^{Cyto} \\cdot n^{Cyto} - k_{4,Vmax} v_{Ro}^{ERm} (u_{AER}^{ER} - u_A^{Cyto}) &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t u_{AER}^{ER} - D_{AER}^{ER} \\nabla^2 u_{AER}^{ER} &= 0 \\qquad \\text{ in } \\Omega^{ER}\\\\\n", + " D_{AER}^{ER} \\nabla u_{AER}^{ER} \\cdot n^{ER} + k_{4,Vmax} v_{Ro}^{ERm} (u_{AER}^{ER} - u_A^{Cyto}) &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t v_{R}^{ERm} - D_{R}^{ERm} \\nabla^2 v_{R}^{ERm} - \n", + " k_{3f} v_R^{ERm} u_B^{Cyto} + k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm}\\\\\n", + " \\nonumber \\\\\n", + " \\partial_t v_{Ro}^{ERm} - D_{Ro}^{ERm} \\nabla^2 v_{Ro}^{ERm} +\n", + " k_{3f} v_R^{ERm} u_B^{Cyto} - k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm}\\\\\n", + "\\end{align}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "28f77cbf", + "metadata": {}, + "source": [ + "## Code imports and initialization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc398816", + "metadata": {}, + "outputs": [], + "source": [ + "import logging\n", + "\n", + "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, mesh, model, mesh_tools, visualization\n", + "from smart.model_assembly import (\n", + " Compartment,\n", + " Parameter,\n", + " Reaction,\n", + " Species,\n", + " SpeciesContainer,\n", + " ParameterContainer,\n", + " CompartmentContainer,\n", + " ReactionContainer,\n", + ")\n", + "from smart.units import unit" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c8650536", + "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": "9ed0899a", + "metadata": {}, + "outputs": [], + "source": [ + "logger = logging.getLogger(\"smart\")\n", + "logger.setLevel(logging.INFO)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "defc1095", + "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": [ + "First, we define the various units for the inputs" + ] + }, + { + "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", + "# Aliases - units used in model\n", + "D_unit = um**2 / sec\n", + "flux_unit = molecule / (um**2 * sec)\n", + "vol_unit = uM\n", + "surf_unit = molecule / um**2" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "46582d26", + "metadata": {}, + "source": [ + "## Generate model\n", + "Next we generate the model described in the equations above.\n", + "\n", + "### Compartments\n", + "As described above, the three compartments are the cytosol (\"Cyto\"), the plasma membrane (\"PM\"), the ER membrane (\"ERm\"), and the ER interior volume (\"ER\").\n", + "\n", + "Note that, as shown, we can also specify nonadjacency for compartments; this is not strictly necessary, but will generally speed up the simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2f34b3c", + "metadata": {}, + "outputs": [], + "source": [ + "Cyto = Compartment(\"Cyto\", 3, um, 1)\n", + "PM = Compartment(\"PM\", 2, um, 10)\n", + "ER = Compartment(\"ER\", 3, um, 2)\n", + "ERm = Compartment(\"ERm\", 2, um, 12)\n", + "PM.specify_nonadjacency(['ERm', 'ER'])\n", + "ERm.specify_nonadjacency(['PM'])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d93ce862", + "metadata": {}, + "source": [ + "Initialize a compartment container and add the 4 compartments to it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "701577e8", + "metadata": {}, + "outputs": [], + "source": [ + "cc = CompartmentContainer()\n", + "cc.add([ERm, ER, PM, Cyto])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "0a6acf0b", + "metadata": {}, + "source": [ + "### Species\n", + "In this case, we have 5 species across 3 different compartments. We define each in turn:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d991278", + "metadata": {}, + "outputs": [], + "source": [ + "A = Species(\"A\", 0.01, vol_unit, 1.0, D_unit, \"Cyto\")\n", + "B = Species(\"B\", 0.0, vol_unit, 1.0, D_unit, \"Cyto\")\n", + "AER = Species(\"AER\", 200.0, vol_unit, 5.0, D_unit, \"ER\")\n", + "# Uniform initial condition of R for simplicity\n", + "R1 = Species(\"R1\", 1.0, surf_unit, 0.02, D_unit, \"ERm\")\n", + "R1o = Species(\"R1o\", 0.0, surf_unit, 0.02, D_unit, \"ERm\")\n", + "\n", + "# Add a species to represent Orai1\n", + "Orai1_open = Species(\"Orai1_open\", 0.0, surf_unit, 0.0, D_unit, \"PM\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b60826cc", + "metadata": {}, + "source": [ + "Create species container and add the 5 species objects to it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e86bebaf", + "metadata": {}, + "outputs": [], + "source": [ + "sc = SpeciesContainer()\n", + "sc.add([R1o, R1, AER, B, A, Orai1_open])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "af900a73", + "metadata": {}, + "source": [ + "### Parameters and Reactions\n", + "\n", + "Parameters and reactions are generally defined together, although the order does not strictly matter. We define them in turn as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2120ac6", + "metadata": {}, + "outputs": [], + "source": [ + "# Degradation of B in the cytosol\n", + "k2f = Parameter(\"k2f\", 10, 1 / sec)\n", + "r2 = Reaction(\n", + " \"r2\", [\"B\"], [], param_map={\"on\": \"k2f\"}, reaction_type=\"mass_action_forward\"\n", + ")\n", + "\n", + "# Activating receptors on ERm with B\n", + "k3f = Parameter(\"k3f\", 100, 1 / (uM * sec))\n", + "k3r = Parameter(\"k3r\", 100, 1 / sec)\n", + "r3 = Reaction(\"r3\", [\"B\", \"R1\"], [\"R1o\"], {\"on\": \"k3f\", \"off\": \"k3r\"})\n", + "# Release of A from ERm to cytosol\n", + "k4Vmax = Parameter(\"k4Vmax\", 2000, 1 / (uM * sec))\n", + "r4 = Reaction(\n", + " \"r4\",\n", + " [\"AER\"],\n", + " [\"A\"],\n", + " param_map={\"Vmax\": \"k4Vmax\"},\n", + " species_map={\"R1o\": \"R1o\", \"uER\": \"AER\", \"u\": \"A\"},\n", + " eqn_f_str=\"Vmax*R1o*(uER-u)\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "7fdd58b8", + "metadata": {}, + "source": [ + "We define an additional reaction as the time-dependent production of species B at the plasma membrane. In this case, we define a pulse-type function as the derivative of an arctan function. Note that this is useful because we can provide an expression to use for pre-integration.\n", + "\n", + "$$\n", + "j_{int}[t] = V_{max} \\arctan\\left({m (t - t_0)}\\right)\\\\\n", + "j_1[t] = \\frac{dj_{int}[t]}{dt} = \\frac{m V_{max}}{1 + m^2 (t-t_0)^2}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51abf270", + "metadata": {}, + "outputs": [], + "source": [ + "Vmax, t0, m = 500, 0.1, 200\n", + "t = sym.symbols(\"t\")\n", + "pulseI = Vmax * sym.atan(m * (t - t0))\n", + "pulse = sym.diff(pulseI, t)\n", + "j1pulse = Parameter.from_expression(\n", + " \"j1pulse\", pulse, flux_unit, use_preintegration=True, preint_sym_expr=pulseI\n", + ")\n", + "r1 = Reaction(\n", + " \"r1\",\n", + " [],\n", + " [\"B\"],\n", + " param_map={\"J\": \"j1pulse\"},\n", + " eqn_f_str=\"J\",\n", + " explicit_restriction_to_domain=\"PM\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b5e218b0", + "metadata": {}, + "source": [ + "We can plot the time-dependent input by converting the sympy expression to a numpy function using lambdify." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0829a6b", + "metadata": {}, + "outputs": [], + "source": [ + "from sympy.utilities.lambdify import lambdify\n", + "pulse_func = lambdify(t, pulse, 'numpy') # returns a numpy-ready function\n", + "tArray = np.linspace(0, 1, 100)\n", + "pulse_vals = pulse_func(tArray)\n", + "plt.plot(tArray, pulse_vals)" + ] + }, + { + "cell_type": "markdown", + "id": "55c27105", + "metadata": {}, + "source": [ + "Finally, we define the store-operated calcium (or here, A) influx. \n", + "This is not immediately straightforward, as a flux at the PM depends on concentration of A in the ER.\n", + "One option would be to explicitly introduce contacts between the ER and PM, but this would imply direct entry of calcium into the ER through Orai1.\n", + "Instead, we assume the likelihood of STIM-Orai1 binding scales with PM-ER distance (if they are very close, membrane fluctuations are likely to allow contact between the ERM and PM occasionally, for instance).\n", + "To adopt this strategy, we must construct a map between each point on the PM and the closest point on the ER, which we do after initializing the model below.\n", + "For now, we define the flux as a parameter varying over the PM mesh by calling `Parameter.mesh_quantity()`\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a56bbccb", + "metadata": {}, + "outputs": [], + "source": [ + "PSOCE = Parameter.mesh_quantity(\"PSOCE\", 0, unit.dimensionless, compartment=\"PM\")\n", + "J0_SOCE = Parameter(\"J0_SOCE\", 1e4, flux_unit)\n", + "\n", + "tauSOCE = Parameter(\"tauSOCE\", 1.0, sec)\n", + "rhoTot = Parameter(\"rhoTot\", 1.0, surf_unit)\n", + "rSOCEOpen = Reaction(\n", + " \"rSOCEOpen\",\n", + " [],\n", + " [\"Orai1_open\"],\n", + " param_map={\"P\": \"PSOCE\", \"rhoTot\": \"rhoTot\", \"tauSOCE\": \"tauSOCE\"},\n", + " eqn_f_str=\"(rhoTot*P - Orai1_open)/tauSOCE\",\n", + " explicit_restriction_to_domain=\"PM\",\n", + ")\n", + "\n", + "rSOCEFlux = Reaction(\n", + " \"rSOCEFlux\",\n", + " [],\n", + " [\"A\"],\n", + " param_map={\"P\": \"PSOCE\", \"J0\": \"J0_SOCE\"},\n", + " eqn_f_str=\"P*J0\",\n", + " explicit_restriction_to_domain=\"PM\",\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "3d883c6e", + "metadata": {}, + "source": [ + "Create containers for parameters and reactions and add all the parameters and reaction objects to them." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ae2c2c1", + "metadata": {}, + "outputs": [], + "source": [ + "pc = ParameterContainer()\n", + "rc = ReactionContainer()\n", + "pc.add([k4Vmax, k3r, k3f, k2f, j1pulse, PSOCE, J0_SOCE, tauSOCE, rhoTot])\n", + "rc.add([r1, r2, r3, r4, rSOCEFlux, rSOCEOpen])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "15c35d39", + "metadata": {}, + "source": [ + "## Create/load in mesh\n", + "\n", + "In SMART we have different levels of meshes:\n", + "- Parent mesh: contains the entire geometry of the problem, including all surfaces and volumes\n", + "- Child meshes: submeshes (sections of the parent mesh) associated with individual compartments. Here, the child meshes are:\n", + " - Cyto: the portion of the outer cube outside of the inner cube, defined by `cell_markers = 1`\n", + " - ER: the inside portion of the inner cube, defined by `cell_markers = 2`\n", + " - PM: surface mesh where x=0, defined by `facet_markers = 10`\n", + " - ERm: surface mesh corresponding to all faces of the inner cube, defined by `facet_markers = 12`\n", + "\n", + "Here we create a UnitCube mesh as the Parent mesh, defined by\n", + "\n", + "$$\n", + "\\Omega = [0, 1] \\times [0, 1] \\times [0, 1] \\subset \\mathbb{R}^3\n", + "$$\n", + "\n", + "\n", + "The ER is a cube within the exterior cube, with dimensions 0.4 by 0.4 by 0.4 and a tunable gap between the ER and PM to test for different strengths of SOCE below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e77f97fe", + "metadata": {}, + "outputs": [], + "source": [ + "ER_PM_gap = 0.05 # must be greater than 0 and less than 0.6\n", + "ER_PM_gap = 0.05 * round(ER_PM_gap/0.05) # round to ensure flat edges in mesh\n", + "def cur_cube_condition(cell, xmin=ER_PM_gap, xmax=ER_PM_gap+0.4):\n", + " \"\"\"\n", + " Returns true when inside an inner cube region defined as:\n", + " xmin <= x <= xmax, ymin <= y <= ymax, zmin <= z <= zmax\n", + " \"\"\"\n", + " ymin = 0.3\n", + " ymax = 0.7\n", + " zmin = 0.3\n", + " zmax = 0.7\n", + " return (\n", + " (xmin - d.DOLFIN_EPS < cell.midpoint().x() < xmax + d.DOLFIN_EPS)\n", + " and (ymin - d.DOLFIN_EPS < cell.midpoint().y() < ymax + d.DOLFIN_EPS)\n", + " and (zmin - d.DOLFIN_EPS < cell.midpoint().z() < zmax + d.DOLFIN_EPS)\n", + " )\n", + "# domain, facet_markers, cell_markers = mesh_tools.create_cubes(condition=cur_cube_condition, N=20)\n", + "domain, facet_markers, cell_markers = mesh_tools.create_nested_substrate(hEdge=5.0)\n", + "visualization.plot_dolfin_mesh(domain, cell_markers, clip_plane=(1,\n", + " 1, 0), clip_origin=(0.5, 0.5, 0.5))\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "035113ed", + "metadata": {}, + "source": [ + "By default, `smart.mesh_tools.create_cubes` marks all faces of the outer cube as \"10\", our marker value associated with PM. Here, since we are only treating the x=0 face as PM, we alter the facet markers on all other faces, setting them equal to zero. They are then treated as no-flux boundaries not belonging to a designated surface compartment." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe56e162", + "metadata": {}, + "outputs": [], + "source": [ + "for face in d.faces(domain):\n", + " if face.midpoint().x() > d.DOLFIN_EPS and facet_markers[face] == 10:\n", + " facet_markers[face] = 0" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c17ebcbb", + "metadata": {}, + "source": [ + "We now save the mesh as an h5 file and then read it into SMART as a `ParentMesh` object. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5657aab1", + "metadata": {}, + "outputs": [], + "source": [ + "mesh_folder = pathlib.Path(\"mesh\")\n", + "mesh_folder.mkdir(exist_ok=True)\n", + "mesh_path = mesh_folder / \"DemoCuboidsMesh.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", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "0f1cf8a4", + "metadata": {}, + "source": [ + "## Model and solver initialization\n", + "\n", + "Now we are ready to set up the model. First we load the default configurations and set some configurations for the current solver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c1a2a92", + "metadata": {}, + "outputs": [], + "source": [ + "conf = config.Config()\n", + "conf.solver.update(\n", + " {\n", + " \"final_t\": 1,\n", + " \"initial_dt\": 0.01,\n", + " \"time_precision\": 6,\n", + " }\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "1511acc7", + "metadata": {}, + "source": [ + "We create a model using the different containers and the parent mesh. For later reference, we save the model information as a pickle file. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3c3c30f", + "metadata": {}, + "outputs": [], + "source": [ + "model_cur = model.Model(pc, sc, cc, rc, conf, parent_mesh)\n", + "model_cur.to_pickle('model_cur.pkl')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "dfb70c74", + "metadata": {}, + "source": [ + "Note that we could later load the model information from the pickle file using the line:\n", + "```Python\n", + "model_cur = model.from_pickle(model_cur.pkl)\n", + "```\n", + "\n", + "Next we need to initialize the model and solver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8976aa2", + "metadata": {}, + "outputs": [], + "source": [ + "model_cur.initialize()" + ] + }, + { + "cell_type": "markdown", + "id": "70f79003", + "metadata": {}, + "source": [ + "Now create mapping between PM and closest point on the ER surface for spatially dependent SOCE." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db8d9dca", + "metadata": {}, + "outputs": [], + "source": [ + "Pfun = PSOCE.dolfin_function # function associated with PSOCE values\n", + "V_PM = Pfun.function_space() # PM function space\n", + "# pull coordinate lists for PM and ERM\n", + "PMcoord = V_PM.tabulate_dof_coordinates()\n", + "ERMcoord = ERm.dolfin_mesh.coordinates()\n", + "# initialize vectors for distances and AER indices\n", + "AER_idx = np.zeros_like(Pfun.vector())\n", + "dist_vals = np.zeros_like(Pfun.vector())\n", + "# get parent mesh and vertex mapping to meshviews\n", + "mesh_ref = model_cur.parent_mesh.dolfin_mesh\n", + "ERM_map = ERm.dolfin_mesh.topology().mapping()[mesh_ref.id()].vertex_map()\n", + "ER_map = ER.dolfin_mesh.topology().mapping()[mesh_ref.id()].vertex_map()\n", + "# find closest ERM point for each PM coordinate\n", + "for i in range(len(PMcoord)):\n", + " curCoord = PMcoord[i]\n", + " dists = np.sqrt((curCoord[0]-ERMcoord[:,0])**2 + \n", + " (curCoord[1]-ERMcoord[:,1])**2 + \n", + " (curCoord[2]-ERMcoord[:,2])**2)\n", + " ERMidx = np.argmin(dists)\n", + " dist_vals[i] = dists[ERMidx]\n", + " global_idx = ERM_map[ERMidx]\n", + " ER_idx = np.nonzero(np.array(ER_map)==global_idx)[0]\n", + " if len(ER_idx) != 1:\n", + " raise ValueError(\"Node not found in ER mesh\")\n", + " AER_idx[i] = d.vertex_to_dof_map(AER.V)[ER_idx][0]\n", + "# now define PSOCE function values according to dist vals and AER values\n", + "d0 = 0.05\n", + "cref = 100.0\n", + "vals_new = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n", + "PSOCE.dolfin_function.vector()[:] = vals_new\n", + "PSOCE.dolfin_function.vector().apply(\"insert\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d05a1a75", + "metadata": {}, + "source": [ + "## Solving system and storing data\n", + "\n", + "We create some XDMF files where we will store the output " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5df9509e", + "metadata": {}, + "outputs": [], + "source": [ + "# Write initial condition(s) to file\n", + "results = dict()\n", + "result_folder = pathlib.Path(f\"results_ER_PM_gap_{ER_PM_gap}\")\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", + "SOCEfile = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"PSOCE.xdmf\"))\n", + "SOCEfile.parameters[\"flush_output\"] = True\n", + "SOCEfile.write(PSOCE.dolfin_function, model_cur.t)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "183376e6", + "metadata": {}, + "source": [ + "We now run the the solver and store the data at each time point to the initialized files. We also integrate A over the cytosolic volume at each time step to monitor the elevation in A over time and we display the concentration of A in the cytosol at t = 0.2 s." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31d42278", + "metadata": {}, + "outputs": [], + "source": [ + "# Set loglevel to warning in order not to pollute notebook output\n", + "logger.setLevel(logging.WARNING)\n", + "\n", + "avg_A = [A.initial_condition]\n", + "# define integration measure and total volume for computing average A at each time point\n", + "dx = d.Measure(\"dx\", domain=model_cur.cc['Cyto'].dolfin_mesh)\n", + "volume = d.assemble_mixed(1.0*dx)\n", + "# Solve\n", + "displayed = False\n", + "while model_cur.t < model_cur.final_t:\n", + " # Solve the system\n", + " model_cur.monolithic_solve()\n", + " # Update PSOCE\n", + " PSOCE.dolfin_function.vector()[:] = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n", + " PSOCE.dolfin_function.vector().apply(\"insert\")\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", + " SOCEfile.write(PSOCE.dolfin_function, model_cur.t)\n", + " # compute average A concentration at each time step\n", + " int_val = d.assemble_mixed(model_cur.sc['A'].u['u']*dx)\n", + " avg_A.append(int_val / volume)\n", + " if model_cur.t >= 0.2 and not displayed:\n", + " visualization.plot(model_cur.sc['A'].u['u'],\n", + " clip_plane=(1, 1, 0), clip_origin=(0.5, 0.5, 0.5))\n", + " displayed = True" + ] + }, + { + "cell_type": "markdown", + "id": "621600d7", + "metadata": {}, + "source": [ + "Finally, we plot the average concentration of A over time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8651e71", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(model_cur.tvec, avg_A, label='PM-ER gap = 50 nm')\n", + "other_case = np.loadtxt('/root/shared/gitrepos/smart-dev/examples/example5/results_ER_PM_gap_0.2/Avals.txt')\n", + "plt.plot(other_case[0,:], other_case[1,:], label='PM-ER gap = 200 nm')\n", + "plt.xlabel('Time (s)')\n", + "plt.ylabel('Cytosolic concentration of A (μM)')\n", + "plt.legend()\n", + "np.savetxt(str(result_folder / f\"Avals.txt\"), [model_cur.tvec, avg_A])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d302bf3f", + "metadata": {}, + "source": [ + "We also compare the AUC for A with previous numerical simulations (regression test)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8b3c866", + "metadata": {}, + "outputs": [], + "source": [ + "tvec = np.zeros(len(model_cur.tvec))\n", + "for i in range(len(model_cur.tvec)):\n", + " tvec[i] = float(model_cur.tvec[i])\n", + "auc_cur = np.trapz(np.array(avg_A), tvec)\n", + "auc_compare = 4.646230684534995\n", + "percent_error = 100*np.abs(auc_cur - auc_compare)/auc_compare\n", + "assert percent_error < .01,\\\n", + " f\"Failed regression test: Example 5 results deviate {percent_error:.3f}% from the previous numerical solution\"" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/example7/example7.ipynb b/examples/example7/example7.ipynb new file mode 100644 index 0000000..a667dd0 --- /dev/null +++ b/examples/example7/example7.ipynb @@ -0,0 +1,336 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "f65f18d7", + "metadata": {}, + "source": [ + "# Example 7: Reaction-diffusion of molecule in a multicellular network\n", + "\n", + "Here, we consider the diffusion of a molecule from a central cell in a multicellular mesh to other cells in the environment." + ] + }, + { + "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\", 3, um, 1)\n", + "PM1 = Compartment(\"PM1\", 2, um, 11) # source cell\n", + "PM2 = Compartment(\"PM2\", 2, um, 12) # other cells\n", + "PM1.specify_nonadjacency(['PM2'])\n", + "PM2.specify_nonadjacency(['PM1'])\n", + "\n", + "cc = CompartmentContainer()\n", + "cc.add([EC, PM1, PM2])\n", + "\n", + "G = Species(\"G\", 1.0, vol_unit, 1000.0, D_unit, \"EC\")\n", + "Rbound = Species(\"Rbound\", 0.0, surf_unit, 0.0, D_unit, \"PM2\")\n", + "sc = SpeciesContainer()\n", + "sc.add([G, Rbound])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "3c56e840", + "metadata": {}, + "source": [ + "Define parameters and reactions, then place in respective containers. Here, there are 3 reactions:\n", + "* r1: release of G from PM1\n", + "* r2: binding of G to PM2\n", + "* r3: degradation of G" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e1f6882", + "metadata": {}, + "outputs": [], + "source": [ + "# G release from cell 1\n", + "j1pulse = Parameter(\"j1pulse\", 1000.0, flux_unit)\n", + "r1 = Reaction(\n", + " \"r1\",\n", + " [],\n", + " [\"G\"],\n", + " param_map={\"J\": \"j1pulse\"},\n", + " eqn_f_str=\"J\",\n", + " explicit_restriction_to_domain=\"PM1\",\n", + ")\n", + "# G binding to PM2\n", + "Rtot = Parameter(\"Rtot\", 100.0, surf_unit)\n", + "kbind = Parameter(\"kbind\", 1.0, 1/(uM*sec))\n", + "kunbind = Parameter(\"kunbind\", 0.01, 1/sec)\n", + "r2 = Reaction(\"r2\", [\"G\"], [\"Rbound\"],\n", + " param_map={\"on\":\"kbind\",\"off\":\"kunbind\",\"Rtot\":\"Rtot\"},\n", + " eqn_f_str=\"on*G*(Rtot-Rbound) - off*Rbound\",\n", + " explicit_restriction_to_domain=\"PM2\")\n", + "\n", + "# G degradation\n", + "kdeg = Parameter(\"kdeg\", 0.01, 1/sec)\n", + "r3 = Reaction(\"r3\", [\"G\"], [], param_map={\"k\":\"kdeg\"},\n", + " eqn_f_str=\"k*G\")\n", + "\n", + "pc = ParameterContainer()\n", + "pc.add([j1pulse,Rtot,kbind,kunbind,kdeg])\n", + "rc = ReactionContainer()\n", + "rc.add([r1,r2,r3])" + ] + }, + { + "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": [ + "domain, facet_markers, cell_markers = mesh_tools.create_multicell(cubeSize=100.0, locVec1=[[0,0,0]],\n", + " locVec2=[[-30,-30,-30], [-30,30,-30], [-30,-30,30], [-30,30,30], \n", + " [30,-30,-30], [30,30,-30], [30,-30,30], [30,30,30]],\n", + " cellRad1 = 10.0, cellRad2 = 10.0, hCube = 5.0, hCell = 2.0)\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 / \"multicell_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", + "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n", + "config_cur.solver.update(\n", + " {\n", + " \"final_t\": 10.0,\n", + " \"initial_dt\": 0.001,\n", + " \"time_precision\": 8,\n", + " \"reset_timestep_for_negative_solution\": True,\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\")\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", + "locVec2 = [[-30,-30,-30], [-30,30,-30], [-30,-30,30], [-30,30,30], \n", + " [30,-30,-30], [30,30,-30], [30,-30,30], [30,30,30]]\n", + "cell_id_fcn = d.MeshFunction(\"size_t\", cc[\"PM2\"].dolfin_mesh, 2, 0)\n", + "for c in d.cells(cc[\"PM2\"].dolfin_mesh):\n", + " xCur = c.midpoint().x()\n", + " yCur = c.midpoint().y()\n", + " zCur = c.midpoint().z()\n", + " for i in range(len(locVec2)):\n", + " print(f\"Checking cell {i+1}, facet {c.index()}\")\n", + " Rcur = np.sqrt((locVec2[i][0]-xCur)**2 + \n", + " (locVec2[i][1]-yCur)**2 +\n", + " (locVec2[i][2]-zCur)**2)\n", + " if Rcur < 11.0:\n", + " cell_id_fcn[c] = i+1\n", + "dx_PM2 = d.Measure(\"dx\", domain=cc[\"PM2\"].dolfin_mesh, subdomain_data=cell_id_fcn)\n", + "avg_concs = []\n", + "avg_conc_cur = []\n", + "sas = []\n", + "for i in range(len(locVec2)):\n", + " sas.append(d.assemble_mixed(1.0*dx_PM2(i+1)))\n", + " avg_conc_cur.append(d.assemble_mixed(sc[\"G\"].sol*dx_PM2(i+1))/sas[-1])\n", + "avg_concs.append(avg_conc_cur)\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", + " avg_conc_cur = []\n", + " for i in range(len(locVec2)):\n", + " avg_conc_cur.append(d.assemble_mixed(sc[\"G\"].sol*dx_PM2(i+1))/sas[i])\n", + " avg_concs.append(avg_conc_cur)\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", + "np.savetxt(\"test_data.txt\", avg_concs)\n", + "for i in range(len(locVec2)):\n", + " plt.plot(model_cur.tvec, np.array(avg_concs)[:,i])\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/examples/example8/example8.ipynb b/examples/example8/example8.ipynb new file mode 100644 index 0000000..e1058fa --- /dev/null +++ b/examples/example8/example8.ipynb @@ -0,0 +1,299 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "f65f18d7", + "metadata": {}, + "source": [ + "# Example 8: Advection-diffusion in a cylindrical geometry\n", + "\n", + "Here, we consider release of a molecule from a cylindrical wall with Pouiselle flow." + ] + }, + { + "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": [ + "Cyto = Compartment(\"Cyto\", 3, um, 1)\n", + "PM = Compartment(\"PM\", 2, um, 10,\n", + " vel=[\"0\",\"0\",\"100.0*[1-(x[0]**2 + x[1]**2//4)]\"])\n", + "\n", + "cc = CompartmentContainer()\n", + "cc.add([Cyto, PM])\n", + "\n", + "A = Species(\"A\", \"1+0.1*z\", vol_unit, 1.0, D_unit, \"Cyto\")\n", + "sc = SpeciesContainer()\n", + "sc.add([A])" + ] + }, + { + "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": [ + "# flux to keep constant \n", + "Awall = Parameter.from_expression(\"Awall\", \"1 + 0.1*z\", vol_unit)\n", + "Jwall = Parameter(\"Jwall\", 1000, flux_unit / vol_unit)\n", + "r1 = Reaction(\n", + " \"r1\",\n", + " [],\n", + " [\"A\"],\n", + " param_map={\"Awall\": \"Awall\", \"Jwall\": \"Jwall\"},\n", + " eqn_f_str=\"Jwall*(Awall-A)\",\n", + " explicit_restriction_to_domain=\"PM\",\n", + ")\n", + "\n", + "pc = ParameterContainer()\n", + "pc.add([Awall, Jwall])\n", + "rc = ReactionContainer()\n", + "rc.add([r1])" + ] + }, + { + "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": [ + "domain, facet_markers, cell_markers = mesh_tools.create_cylinders(outerRad = 2.0, innerRad=0, outerLength=10.0, innerLength=0,\n", + " hEdge=0.5)\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", + "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n", + "config_cur.solver.update(\n", + " {\n", + " \"final_t\": 100.0,\n", + " \"initial_dt\": 0.01,\n", + " \"time_precision\": 8,\n", + " \"reset_timestep_for_negative_solution\": True,\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\")\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,)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b755f19b", + "metadata": {}, + "outputs": [], + "source": [ + "u = sc[\"A\"].sol\n", + "udiff = u - d.Expression(\"1.0+0.1*x[2]\")" + ] + } + ], + "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/pyproject.toml b/pyproject.toml index 0bbb4e3..4a5e39a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ requires = ["setuptools>=61.0.0", "wheel"] [project] name = "fenics-smart" -version = "2.2.3" +version = "2.3.1.alpha.1" description = "Spatial Modeling Algorithms for Reactions and Transport (SMART) is a high-performance finite-element-based simulation package for model specification and numerical simulation of spatially-varying reaction-transport processes in biological cells." authors = [{ name = "Justin Laughlin", email = "justinglaughlin@gmail.com" }] license = { file = "LICENSE" } diff --git a/smart/mesh.py b/smart/mesh.py index 96774ce..ba93abd 100644 --- a/smart/mesh.py +++ b/smart/mesh.py @@ -152,6 +152,10 @@ class ParentMesh(_Mesh): as an xml or xdmf file containing a vertex mesh function of type double extra_keys (list): list of names of extra keys to load from hdf5 mesh file, specifying other subdomains besides main compartments (optional) + key_for_cells (string): a string specifying the extra key to use as a + mesh function to define cell markers in mesh (optional) + key_for_facets (string): a string specifying the extra key to use as a + mesh function to define facet markers in mesh (optional) """ mesh_filename: str @@ -161,6 +165,8 @@ class ParentMesh(_Mesh): use_partition: bool curvature: d.MeshFunction extra_keys: list = [] + key_for_cells: str = "" + key_for_facets: str = "" def __init__( self, @@ -171,6 +177,8 @@ def __init__( mpi_comm=d.MPI.comm_world, curvature=None, extra_keys=[], + key_for_cells="", + key_for_facets="", ): super().__init__(name) self.use_partition = use_partition @@ -202,6 +210,30 @@ def __init__( # Otherwise just take what we got self.curvature = curvature + # set cells and/or facets according to extra keys if applicable + if key_for_cells != "": + try: + idx = self.extra_keys.index(key_for_cells) + except ValueError: + raise ValueError(f"'{key_for_cells}' does not match an extra key") + assert self.subdomains[idx].dim() == self.dimensionality, ( + f"Mesh function associated with '{key_for_cells}' " + "does not match mesh cell dimension" + ) + self.mf["cells"] = self.subdomains[idx] + logger.info(f"Cell mesh function loaded from key '{key_for_cells}'") + if key_for_facets != "": + try: + idx = self.extra_keys.index(key_for_facets) + except ValueError: + raise ValueError(f"'{key_for_cells}' does not match an extra key") + assert self.subdomains[idx].dim() == self.dimensionality - 1, ( + f"Mesh function associated with '{key_for_facets}' " + "does not match mesh facet dimension" + ) + self.mf["facets"] = self.subdomains[idx] + logger.info(f"Facet mesh function loaded from key '{key_for_facets}'") + def get_mesh_from_id(self, id): "Find the mesh that has the matching id." # find the mesh in that has the matching id @@ -290,8 +322,9 @@ def read_parent_mesh_functions_from_file(self): assert len(self.child_meshes) > 0 # Init mesh functions - self.mf["cells"] = self._read_parent_mesh_function_from_file(volume_dim) - if self.has_surface: + if "cells" not in self.mf.keys(): + self.mf["cells"] = self._read_parent_mesh_function_from_file(volume_dim) + if self.has_surface and "facets" not in self.mf.keys(): self.mf["facets"] = self._read_parent_mesh_function_from_file(surface_dim) # If any cell markers are given as a list we also create mesh diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py index e1a700a..63fe61b 100644 --- a/smart/mesh_tools.py +++ b/smart/mesh_tools.py @@ -719,6 +719,254 @@ def meshSizeCallback(dim, tag, x, y, z, lc): return (dmesh, mf2, mf3) +def create_multicell( + cubeSize: float = 100.0, + locVec1: list = [[0, 0, 0]], + locVec2: list = [], + cellRad1: float = 10.0, + cellRad2: float = 10.0, + hCube: float = 0, + hCell: float = 0, + interface_marker1: int = 11, + interface_marker2: int = 12, + outer_marker: int = 10, + extracell_tag: int = 1, + cell_vol_tag1: int = 2, + cell_vol_tag2: int = 3, + comm: MPI.Comm = d.MPI.comm_world, + verbose: bool = False, +) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]: + """ + Creates a mesh with an outer cube containing embedded cells at specified locations. + Cells can be of two types, type 1 and type 2 throughout. + Args: + cubeSize: Length of cube sides + locVec1: list of cell 1 locations (each list element is a list [x,y,z]) + locVec2: list of cell 2 locations (each list element is a list [x,y,z]) + cellRad1: either a float, a list of floats, or a list of lists, each + containing a list of three floats for each axis [rx,ry,rz] + cellRad2: either a float, a list of floats, or a list of lists, each + containing a list of three floats for each axis [rx,ry,rz] + hCube: maximum mesh size for cube + hCell: maximum mesh size for cell surface + interface_marker1: The value to mark facets on cell 1 surface + outer_marker: The value to mark facets on the surface of the cube + extracell_tag: Tag for extracellular volume + cell_vol_tag1: Tag for the volume of type 1 cells + outer_vol_tag2: Tag for the volume of type 2 cells + comm: MPI communicator to create the mesh with + verbose: If true print gmsh output, else skip + Returns: + A triplet (mesh, facet_marker, cell_marker) + """ + import gmsh + + if np.isclose(cubeSize, 0): + ValueError("Outer cube size is equal to zero") + if np.isclose(hCube, 0): + hCube = 0.1 * max(cubeSize) + + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", int(verbose)) + + gmsh.model.add("multicell") + # first add outer cube + cube = gmsh.model.occ.addBox( + -cubeSize / 2, -cubeSize / 2, -cubeSize / 2, cubeSize, cubeSize, cubeSize + ) + meanRad1 = 0 + meanRad2 = 0 + if (len(locVec1) == 0) and (len(locVec2) == 0): + # Just a cube! + gmsh.model.occ.synchronize() + gmsh.model.add_physical_group(3, [cube], tag=extracell_tag) + facets = gmsh.model.getBoundary([(3, cube)]) + gmsh.model.add_physical_group(2, [facets[0][1]], tag=outer_marker) + else: + # Add cells + cell_list = [] + cellRad1Vec = [] + cellRad2Vec = [] + # first add source cell(s) + for i in range(len(locVec1)): + if isinstance(cellRad1, float) or isinstance(cellRad1, int): + cellRadCur = float(cellRad1) + elif len(locVec1) == 1 and len(cellRad1) == 3: + cellRadCur = cellRad1 + elif len(cellRad1) == len(locVec1): + cellRadCur = cellRad1[i] + else: + raise ValueError("Radii must be floats or lists of 3 floats") + if isinstance(cellRadCur, float): + cellRads = [cellRadCur, cellRadCur, cellRadCur] + if len(cellRadCur) == 3: + cellRads = cellRadCur + else: + raise ValueError("Radii must be floats or lists of 3 floats") + + cur_tag = gmsh.model.occ.addSphere(locVec1[i][0], locVec1[i][1], locVec1[i][2], 1.0) + gmsh.model.occ.dilate( + [(3, cur_tag)], + locVec1[i][0], + locVec1[i][1], + locVec1[i][2], + cellRads[0], + cellRads[1], + cellRads[2], + ) + cell_list.append((3, cur_tag)) + meanRad1 += np.mean(cellRads) + cellRad1Vec.append(cellRads) + if meanRad1 > 0: + meanRad1 /= len(locVec1) + # now add additional cells + for i in range(len(locVec2)): + if isinstance(cellRad2, float) or isinstance(cellRad2, int): + cellRadCur = float(cellRad2) + elif len(locVec2) == 1 and len(cellRad2) == 3: + cellRadCur = cellRad2 + elif len(cellRad2) == len(locVec2): + cellRadCur = cellRad2[i] + else: + raise ValueError("Radii must be floats or lists of 3 floats") + if isinstance(cellRadCur, float): + cellRads = [cellRadCur, cellRadCur, cellRadCur] + if len(cellRadCur) == 3: + cellRads = cellRadCur + else: + raise ValueError("Radii must be floats or lists of 3 floats") + + cur_tag = gmsh.model.occ.addSphere(locVec2[i][0], locVec2[i][1], locVec2[i][2], 1.0) + gmsh.model.occ.dilate( + [(3, cur_tag)], + locVec2[i][0], + locVec2[i][1], + locVec2[i][2], + cellRads[0], + cellRads[1], + cellRads[2], + ) + cell_list.append((3, cur_tag)) + meanRad2 = np.mean(cellRads) + cellRad2Vec.append(cellRads) + if meanRad2 > 0: + meanRad2 /= len(locVec2) + # define hCell if it was previously set to zero + if np.isclose(hCell, 0): + hCell = ( + 0.2 * cubeSize + if (meanRad1 == 0 and meanRad2 == 0) + else 0.2 * max([meanRad1, meanRad2]) + ) + # Create interface between cells and extracell + full_geo, maps = gmsh.model.occ.fragment([(3, cube)], cell_list) + cube_map = maps[0] + cell_maps = maps[1:] + gmsh.model.occ.synchronize() + + # Get the outer boundary + outer_shells = gmsh.model.getBoundary(full_geo, oriented=False) + # Get the inner boundary + inner_shells1 = [] + inner_shells2 = [] + for i in range(0, len(locVec1)): + inner_shells1.append(gmsh.model.getBoundary(cell_maps[i], oriented=False)) + for i in range(len(locVec1), len(cell_maps)): + inner_shells2.append(gmsh.model.getBoundary(cell_maps[i], oriented=False)) + # Add physical markers for facets + gmsh.model.add_physical_group(2, [faces[1] for faces in outer_shells], tag=outer_marker) + gmsh.model.add_physical_group( + 2, [faces[0][1] for faces in inner_shells1], tag=interface_marker1 + ) + gmsh.model.add_physical_group( + 2, [faces[0][1] for faces in inner_shells2], tag=interface_marker2 + ) + + # Physical markers for all volumes + all_volumes = [tag[1] for tag in cube_map] + inner_volumes1 = [tag[0][1] for tag in cell_maps[0 : len(locVec1)]] + inner_volumes2 = [tag[0][1] for tag in cell_maps[len(locVec1) :]] + outer_volume = [] + for vol in all_volumes: + if (vol not in inner_volumes1) and (vol not in inner_volumes2): + outer_volume.append(vol) + gmsh.model.add_physical_group(3, outer_volume, tag=extracell_tag) + gmsh.model.add_physical_group(3, inner_volumes1, tag=cell_vol_tag1) + gmsh.model.add_physical_group(3, inner_volumes2, tag=cell_vol_tag2) + + def meshSizeCallback(dim, tag, x, y, z, lc): + # mesh length is hEdge at the PM (defaults to 0.1*outerRad, + # or set when calling function) and hInnerEdge at the ERM + # (defaults to 0.2*innerRad, or set when calling function) + # between these, the value is interpolated based on r (polar coord), + # and inside the value is interpolated between hInnerEdge and 0.2*innerRad + # if innerRad=0, then the mesh length is interpolated between + # hEdge at the PM and 0.2*outerRad in the center + + if len(locVec1) == 0 and len(locVec2) == 0: + return hCube + elif len(locVec1) == 0: + cell_locs1 = [np.inf] + cell_locs2 = np.sqrt( + ((x - np.array(locVec2)[:, 0]) / np.array(cellRad2Vec)[:, 0]) ** 2 + + ((y - np.array(locVec2)[:, 1]) / np.array(cellRad2Vec)[:, 1]) ** 2 + + ((z - np.array(locVec2)[:, 2]) / np.array(cellRad2Vec)[:, 2]) ** 2 + ) + elif len(locVec2) == 0: + cell_locs1 = np.sqrt( + ((x - np.array(locVec1)[:, 0]) / np.array(cellRad1Vec)[:, 0]) ** 2 + + ((y - np.array(locVec1)[:, 1]) / np.array(cellRad1Vec)[:, 1]) ** 2 + + ((z - np.array(locVec1)[:, 2]) / np.array(cellRad1Vec)[:, 2]) ** 2 + ) + cell_locs2 = [np.inf] + else: + cell_locs1 = np.sqrt( + ((x - np.array(locVec1)[:, 0]) / np.array(cellRad1Vec)[:, 0]) ** 2 + + ((y - np.array(locVec1)[:, 1]) / np.array(cellRad1Vec)[:, 1]) ** 2 + + ((z - np.array(locVec1)[:, 2]) / np.array(cellRad1Vec)[:, 2]) ** 2 + ) + cell_locs2 = np.sqrt( + ((x - np.array(locVec2)[:, 0]) / np.array(cellRad2Vec)[:, 0]) ** 2 + + ((y - np.array(locVec2)[:, 1]) / np.array(cellRad2Vec)[:, 1]) ** 2 + + ((z - np.array(locVec2)[:, 2]) / np.array(cellRad2Vec)[:, 2]) ** 2 + ) + closest_cell1 = min(cell_locs1) + closest_cell2 = min(cell_locs2) + if (closest_cell1 < 1.0) or (closest_cell2 < 1.0): + return hCell + else: + if closest_cell1 < closest_cell2: + cellWeight = np.exp(-(closest_cell1 - 1.0) / 0.1) + else: + cellWeight = np.exp(-(closest_cell2 - 1.0) / 0.1) + return hCell * cellWeight + hCube * (1 - cellWeight) + + 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(3) + rank = MPI.COMM_WORLD.rank + tmp_folder = pathlib.Path(f"tmp_extracell_{cubeSize}_{cellRad1}_{cellRad2}_{rank}") + tmp_folder.mkdir(exist_ok=True) + gmsh_file = tmp_folder / "extracell.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, 3, comm) + # remove tmp mesh and tmp folder + gmsh_file.unlink(missing_ok=False) + tmp_folder.rmdir() + # return dolfin mesh, mf2 (2d tags) and mf3 (3d tags) + return (dmesh, mf2, mf3) + + def create_ellipses( xrad_outer: float = 3.0, yrad_outer: float = 1.0, @@ -1300,6 +1548,291 @@ def meshSizeCallback(dim, tag, x, y, z, lc): return (dmesh, mf2, mf3) +def create_2Dcell_xy( + outerExpr: str = "", + innerExpr: str = "", + hEdge: float = 0, + hInnerEdge: float = 0, + interface_marker: int = 12, + outer_marker: int = 10, + inner_tag: int = 2, + outer_tag: int = 1, + comm: MPI.Comm = d.MPI.comm_world, + verbose: bool = False, + half_cell: bool = True, + return_curvature: bool = False, +) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]: + """ + Creates a 2D mesh of a cell profile, with the bounding curve defined in + terms of r and z (e.g. unit circle would be "r**2 + (z-1)**2 - 1) + It is assumed that substrate is present at z = 0, so if the curve extends + below z = 0 , there is a sharp cutoff. + If half_cell = True, only have of the contour is constructed, with a + left zero-flux boundary at r = 0. + Can include one compartment inside another compartment. + Recommended for use with the axisymmetric feature of SMART. + + Args: + outerExpr: String implicitly defining an r-z curve for the outer surface + innerExpr: String implicitly defining an r-z curve for the inner surface + 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 edge of the outer ellipse with + inner_tag: The value to mark the inner ellipse surface with + outer_tag: The value to mark the outer ellipse surface with + comm: MPI communicator to create the mesh with + verbose: If true print gmsh output, else skip + half_cell: If true, consider r=0 the symmetry axis for an axisymm shape + Returns: + Tuple (mesh, facet_marker, cell_marker) + """ + import gmsh + + if outerExpr == "": + ValueError("Outer surface is not defined") + + if return_curvature: + # create full mesh for curvature analysis and then map onto half mesh + # if half_cell_with_curvature is True + half_cell_with_curvature = half_cell + half_cell = False + + rValsOuter, zValsOuter = implicit_curve(outerExpr) + + if not innerExpr == "": + rValsInner, zValsInner = implicit_curve(innerExpr) + zMid = np.mean(zValsInner) + ROuterVec = np.sqrt(rValsOuter**2 + (zValsOuter - zMid) ** 2) + RInnerVec = np.sqrt(rValsInner**2 + (zValsInner - zMid) ** 2) + maxOuterDim = max(ROuterVec) + maxInnerDim = max(RInnerVec) + else: + zMid = np.mean(zValsOuter) + ROuterVec = np.sqrt(rValsOuter**2 + (zValsOuter - zMid) ** 2) + maxOuterDim = max(ROuterVec) + if np.isclose(hEdge, 0): + hEdge = 0.1 * maxOuterDim + if np.isclose(hInnerEdge, 0): + hInnerEdge = 0.2 * maxOuterDim if innerExpr == "" else 0.2 * maxInnerDim + # Create the 2D mesh using gmsh + gmsh.initialize() + gmsh.option.setNumber("General.Terminal", int(verbose)) + gmsh.model.add("2DCell") + # first add outer body and revolve + outer_tag_list = [] + for i in range(len(rValsOuter)): + cur_tag = gmsh.model.occ.add_point(rValsOuter[i], zValsOuter[i], 0.0) + outer_tag_list.append(cur_tag) + outer_spline_tag = gmsh.model.occ.add_spline(outer_tag_list) + if not half_cell: + outer_tag_list2 = [] + for i in range(len(rValsOuter)): + cur_tag = gmsh.model.occ.add_point(-rValsOuter[i], zValsOuter[i], 0.0) + outer_tag_list2.append(cur_tag) + outer_spline_tag2 = gmsh.model.occ.add_spline(outer_tag_list2) + if np.isclose(zValsOuter[-1], 0): # then include substrate at z=0 + if half_cell: + origin_tag = gmsh.model.occ.add_point(0, 0, 0) + symm_axis_tag = gmsh.model.occ.add_line(origin_tag, outer_tag_list[0]) + bottom_tag = gmsh.model.occ.add_line(origin_tag, outer_tag_list[-1]) + outer_loop_tag = gmsh.model.occ.add_curve_loop( + [outer_spline_tag, bottom_tag, symm_axis_tag] + ) + else: + bottom_tag = gmsh.model.occ.add_line(outer_tag_list[-1], outer_tag_list2[-1]) + outer_loop_tag = gmsh.model.occ.add_curve_loop( + [outer_spline_tag, outer_spline_tag2, bottom_tag] + ) + else: + if half_cell: + symm_axis_tag = gmsh.model.occ.add_line(outer_tag_list[0], outer_tag_list[-1]) + outer_loop_tag = gmsh.model.occ.add_curve_loop([outer_spline_tag, symm_axis_tag]) + else: + outer_loop_tag = gmsh.model.occ.add_curve_loop([outer_spline_tag, outer_spline_tag2]) + cell_plane_tag = gmsh.model.occ.add_plane_surface([outer_loop_tag]) + + if innerExpr == "": + # No inner shape in this case + gmsh.model.occ.synchronize() + gmsh.model.add_physical_group(2, [cell_plane_tag], tag=outer_tag) + facets = gmsh.model.getBoundary([(2, cell_plane_tag)]) + facet_tag_list = [] + for i in range(len(facets)): + facet_tag_list.append(facets[i][1]) + if half_cell: # if half, set symmetry axis to 0 (no flux) + xmin, ymin, zmin = (-hInnerEdge / 10, -1, -hInnerEdge / 10) + xmax, ymax, zmax = (hInnerEdge / 10, max(zValsOuter) + 1, hInnerEdge / 10) + 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]) + gmsh.model.add_physical_group(1, symm_bound_markers, tag=0) + gmsh.model.add_physical_group(1, facet_tag_list, tag=outer_marker) + else: + # Add inner shape + inner_tag_list = [] + for i in range(len(rValsInner)): + cur_tag = gmsh.model.occ.add_point(rValsInner[i], zValsInner[i], 0.0) + inner_tag_list.append(cur_tag) + inner_spline_tag = gmsh.model.occ.add_spline(inner_tag_list) + if half_cell: + symm_inner_tag = gmsh.model.occ.add_line(inner_tag_list[0], inner_tag_list[-1]) + inner_loop_tag = gmsh.model.occ.add_curve_loop([inner_spline_tag, symm_inner_tag]) + else: + inner_tag_list2 = [] + for i in range(len(rValsInner)): + cur_tag = gmsh.model.occ.add_point(-rValsInner[i], zValsInner[i], 0.0) + inner_tag_list2.append(cur_tag) + inner_spline_tag2 = gmsh.model.occ.add_spline(inner_tag_list2) + inner_loop_tag = gmsh.model.occ.add_curve_loop([inner_spline_tag, inner_spline_tag2]) + 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 + if half_cell: # if half, set symmetry axis to 0 (no flux) + xmin, ymin, zmin = (-hInnerEdge / 10, -1, -hInnerEdge / 10) + xmax, ymax, zmax = (hInnerEdge / 10, max(zValsOuter) + 1, hInnerEdge / 10) + 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_tag) + gmsh.model.add_physical_group(2, inner_volume, tag=inner_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 + rCur = np.sqrt(x**2 + z**2) + RCur = np.sqrt(rCur**2 + (y - zMid) ** 2) + outer_dist = np.sqrt((rCur - rValsOuter) ** 2 + (z - zValsOuter) ** 2) + np.append(outer_dist, z) # include the distance from the substrate + dist_to_outer = min(outer_dist) + if innerExpr == "": + lc3 = 0.2 * maxOuterDim + dist_to_inner = RCur + in_outer = True + else: + inner_dist = np.sqrt((rCur - rValsInner) ** 2 + (y - zValsInner) ** 2) + dist_to_inner = min(inner_dist) + inner_idx = np.argmin(inner_dist) + inner_rad = RInnerVec[inner_idx] + R_rel_inner = RCur / inner_rad + lc3 = max(hInnerEdge, 0.2 * maxInnerDim) + in_outer = R_rel_inner > 1 + lc1 = hEdge + lc2 = hInnerEdge + if in_outer: + lcTest = lc1 + (lc2 - lc1) * (dist_to_outer) / (dist_to_inner + dist_to_outer) + 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, mf1, mf2 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 2, comm) + # remove tmp mesh and tmp folder + gmsh_file.unlink(missing_ok=False) + tmp_folder.rmdir() + # return dolfin mesh, mf1 (1d tags) and mf2 (2d tags) + if return_curvature: + if innerExpr == "": + facet_list = [outer_marker] + cell_list = [outer_tag] + else: + facet_list = [outer_marker, interface_marker] + cell_list = [outer_tag, inner_tag] + if half_cell_with_curvature: # will likely not work in parallel... + dmesh_half, mf1_half, mf2_half = create_2Dcell( + outerExpr, + innerExpr, + hEdge, + hInnerEdge, + interface_marker, + outer_marker, + inner_tag, + outer_tag, + comm, + verbose, + half_cell=True, + return_curvature=False, + ) + kappa_mf = compute_curvature( + dmesh, mf1, mf2, facet_list, cell_list, half_mesh_data=(dmesh_half, mf1_half) + ) + (dmesh, mf1, mf2) = (dmesh_half, mf1_half, mf2_half) + else: + kappa_mf = compute_curvature(dmesh, mf1, mf2, facet_list, cell_list) + return (dmesh, mf1, mf2, kappa_mf) + else: + return (dmesh, mf1, mf2) + + def gmsh_to_dolfin( gmsh_file_name: str, tmp_folder: pathlib.Path = pathlib.Path("tmp_folder"), @@ -1340,6 +1873,8 @@ def gmsh_to_dolfin( # convert cell mesh cells = mesh_in.get_cells_type(cell_type) cell_data = mesh_in.get_cell_data("gmsh:physical", cell_type) # extract values of tags + if dimension == 2 and np.all(mesh_in.points[:, 2] == 0): + mesh_in.points = mesh_in.points[:, :2] # then prune z values out_mesh_cell = meshio.Mesh( points=mesh_in.points, cells={cell_type: cells}, diff --git a/smart/model.py b/smart/model.py index 8f51475..82c2c9a 100644 --- a/smart/model.py +++ b/smart/model.py @@ -756,14 +756,14 @@ def _init_4_0_initialize_dolfin_parameters(self): if parameter.type == ParameterType.constant: parameter.dolfin_constant = d.Constant(parameter.value, name=parameter.name) elif parameter.type == ParameterType.expression and parameter.is_space_dependent: + c_code = sym.printing.ccode(parameter.sym_expr) + c_code = c_code.replace("log(", "std::log(") # use higher degree to avoid interpolation error - parameter.dolfin_expression = d.Expression( - sym.printing.ccode(parameter.sym_expr), t=self.T, degree=3 - ) + parameter.dolfin_expression = d.Expression(c_code, t=self.T, degree=3) elif parameter.type == ParameterType.expression and not parameter.use_preintegration: - parameter.dolfin_expression = d.Expression( - sym.printing.ccode(parameter.sym_expr), t=self.T, degree=1 - ) + c_code = sym.printing.ccode(parameter.sym_expr) + c_code = c_code.replace("log(", "std::log(") + parameter.dolfin_expression = d.Expression(c_code, t=self.T, degree=1) elif parameter.type == ParameterType.expression and parameter.use_preintegration: parameter.dolfin_constant = d.Constant(parameter.value, name=parameter.name) elif parameter.type == ParameterType.from_file: @@ -1035,6 +1035,107 @@ def _init_4_7_set_initial_conditions(self): new_vals = vec_new[mesh_map] # reorder to match dof ordering vec.set_local(new_vals) vec.apply("insert") + elif parameter.type == ParameterType.mesh_quantity: + if parameter.compartment not in self.cc.keys: + raise ValueError( + f"Compartment name {parameter.compartment} for parameter" + f"{parameter.name} does not match a known compartment" + ) + V_cur = d.FunctionSpace(self.cc[parameter.compartment].dolfin_mesh, "P", 1) + parameter.dolfin_function = d.Function(V_cur) + parameter.dolfin_function.vector()[:] = parameter.value + parameter.dolfin_function.vector().apply("insert") + + for compartment in self._active_compartments: + # if vel is string, generate expr for advection + compartment.vel_logic = np.any([vel != 0.0 for vel in compartment.vel]) + compartment.deform_logic = np.any([deform != 0.0 for deform in compartment.deform]) + if compartment.vel_logic and compartment.deform_logic: + raise ValueError("Cannot prescribe both velocity and deformation") + elif compartment.vel_logic: + x, y, z = (sym.Symbol(f"x[{i}]") for i in range(3)) + vel_expr = [None] * len(compartment.vel) + # then this needs to have free symbols inserted + for i in range(len(compartment.vel)): + if isinstance(compartment.vel[i], float): + vel_expr[i] = f"{compartment.vel[i]}" + elif isinstance(compartment.vel[i], str): + vel_str = compartment.vel[i] + # Parse the given string to create a sympy expression + sym_expr = parse_expr(vel_str).subs({"x": x, "y": y, "z": z}) + free_symbols = [str(x) for x in sym_expr.free_symbols] + if not {"x[0]", "x[1]", "x[2]", "t"}.issuperset(free_symbols): + # could add other keywords for spatial dependence in the future + raise NotImplementedError + else: + vel_expr[i] = sym.printing.ccode(sym_expr) + else: + raise NotImplementedError("Velocity must be float or string") + compartment.vel_expr = d.Expression(vel_expr, degree=1, t=self.T) + V_vector = d.VectorFunctionSpace(compartment.dolfin_mesh, "P", 1) + compartment.vel_func = d.interpolate(compartment.vel_expr, V_vector) + elif compartment.deform_logic: + x, y, z = (sym.Symbol(f"x[{i}]") for i in range(3)) + deform_expr = [None] * len(compartment.deform) + # then this needs to have free symbols inserted + for i in range(len(compartment.deform)): + if isinstance(compartment.deform[i], float): + deform_expr[i] = f"{compartment.deform[i]}" + elif isinstance(compartment.deform[i], str): + deform_str = compartment.deform[i] + # Parse the given string to create a sympy expression + sym_expr = parse_expr(deform_str).subs({"x": x, "y": y, "z": z}) + free_symbols = [str(x) for x in sym_expr.free_symbols] + if not {"x[0]", "x[1]", "x[2]", "t"}.issuperset(free_symbols): + # could add other keywords for spatial dependence in the future + raise NotImplementedError + else: + deform_expr[i] = sym.printing.ccode(sym_expr) + else: + raise NotImplementedError("Deformation must be float or string") + compartment.deform_expr = d.Expression(deform_expr, degree=1, t=self.T) + V_vector = d.VectorFunctionSpace(compartment.dolfin_mesh, "P", 1) + compartment.deform_func = d.interpolate(compartment.deform_expr, V_vector) + compartment.deform_prev = d.interpolate(compartment.deform_expr, V_vector) + if not compartment.is_volume: # then is a surface and must compute normals + surf_coords = compartment.dolfin_mesh.coordinates() + for c in self.cc: # find an adjacent volume! + if c.is_volume: + mesh_test = c.dolfin_mesh + test_coords = mesh_test.coordinates() + adjacent = True + for i in range(len(surf_coords)): + if surf_coords[i] not in test_coords: + adjacent = False + break + if adjacent: + mesh_ref = mesh_test + break + if not adjacent: + raise ValueError("Could not find an adjacent volume to compute normals!") + ref_normals = d.FacetNormal(mesh_ref) + Vcur = d.VectorFunctionSpace(mesh_ref, "P", 1) + ucur = d.TrialFunction(Vcur) + vcur = d.TestFunction(Vcur) + ds = d.Measure("ds", mesh_ref) + a = d.inner(ucur, vcur) * ds + lform = d.inner(ref_normals, vcur) * ds + A = d.assemble(a, keep_diagonal=True) + L = d.assemble(lform) + A.ident_zeros() + nh = d.Function(Vcur) + d.solve(A, nh.vector(), L) # project facet normals onto CG1 + Vbound = d.VectorFunctionSpace(compartment.mesh.dolfin_mesh, "P", 1) + norm_calc = d.interpolate(nh, Vbound) + nvec = norm_calc.vector()[:] + # normalize magnitudes for consistency (unit normal) + for i in range(int(len(nvec) / 3)): + vec_cur = nvec[3 * i : 3 * (i + 1)] + mag_cur = np.sqrt(vec_cur[0] ** 2 + vec_cur[1] ** 2 + vec_cur[2] ** 2) + nvec[3 * i : 3 * (i + 1)] = vec_cur / mag_cur + norm_calc.vector().set_local(nvec) + norm_calc.vector().apply("insert") + compartment.normals = norm_calc def _init_5_1_reactions_to_fluxes(self): """Convert reactions to flux objects""" @@ -1102,6 +1203,16 @@ def _init_5_2_create_variational_forms(self): / unit.s * species.compartment.compartment_units**species.compartment.dimensionality ) + lagrange = species.compartment.deform_logic + # define current jacobian + if lagrange: + udef = species.compartment.deform_func + F = d.Identity(3) + d.grad(udef) + J = d.det(F) + else: + J = d.Constant(1.0) + if self.config.flags["axisymmetric_model"]: + J = x[0] * J # diffusion term if species.D == 0: logger.debug( @@ -1110,30 +1221,74 @@ def _init_5_2_create_variational_forms(self): extra=dict(format_type="log"), ) else: - 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 + if lagrange: # then nonzero advection + # alt Lagrangian form + D_constant = d.Constant(D, name=f"D_{species.name}") + Dform = ( + D_constant + * J + * d.inner(d.dot(d.inv(F.T), d.grad(u)), d.dot(d.inv(F.T), d.grad(v))) + * dx + ) + self.forms.add( + Form( + f"diffusion_{species.name}", + Dform, + species, + "diffusion", + Dform_units, + True, + linear_wrt_comp, + ) + ) else: - Dform = D_constant * d.inner(d.grad(u), d.grad(v)) * dx - # exponent is -2 because of two gradients - + D_constant = d.Constant(D, name=f"D_{species.name}") + Dform = J * D_constant * d.inner(d.grad(u), d.grad(v)) * dx + # exponent is -2 because of two gradients + + self.forms.add( + Form( + f"diffusion_{species.name}", + Dform, + species, + "diffusion", + Dform_units, + True, + linear_wrt_comp, + ) + ) + if lagrange: # account for volume changes due to advection + vel_approx = (udef - species.compartment.deform_prev) / self.dT + Aform = J * u * v * d.inner(d.inv(F.T), d.grad(vel_approx)) * dx + # Aform = J * u * v * d.inner(vel, d.dot(d.inv(F.T), d.grad(u))) * dx self.forms.add( Form( - f"diffusion_{species.name}", - Dform, + f"advection_{species.name}", + Aform, species, - "diffusion", - Dform_units, + "advection", + mass_form_units, + True, + linear_wrt_comp, + ) + ) + elif species.compartment.vel_logic: # then nonzero advection + vel = species.compartment.vel_func + Aform = J * u * v * d.div(vel) * dx + J * d.inner(v * vel, d.grad(u)) * dx + self.forms.add( + Form( + f"advection_{species.name}", + Aform, + species, + "advection", + mass_form_units, True, linear_wrt_comp, ) ) # mass (time derivative) terms - if self.config.flags["axisymmetric_model"]: - Muform = x[0] * (u) * v / self.dT * dx - else: - Muform = (u) * v / self.dT * dx + Muform = J * u * v / self.dT * dx self.forms.add( Form( f"mass_u_{species.name}", @@ -1145,10 +1300,7 @@ def _init_5_2_create_variational_forms(self): linear_wrt_comp, ) ) - if self.config.flags["axisymmetric_model"]: - Munform = x[0] * (-un) * v / self.dT * dx - else: - Munform = (-un) * v / self.dT * dx + Munform = J * (-un) * v / self.dT * dx self.forms.add( Form( f"mass_un_{species.name}", @@ -1855,6 +2007,16 @@ def update_time_dependent_parameters(self): dt = float(self.dt) tn = float(self.tn) + # Update velocity or deformation if applicable + for compartment in self._active_compartments: + if compartment.vel_logic and not compartment.manual_update: + Vcur = compartment.vel_func.function_space() + compartment.vel_func.assign(d.interpolate(compartment.vel_expr, Vcur)) + elif compartment.deform_logic and not compartment.manual_update: + Vcur = compartment.deform_func.function_space() + compartment.deform_prev.assign(compartment.deform_func.copy(deepcopy=True)) + compartment.deform_func.assign(d.interpolate(compartment.deform_expr, Vcur)) + # Update time dependent parameters for parameter_name, parameter in self.pc.items: new_value = None @@ -1909,9 +2071,9 @@ def update_time_dependent_parameters(self): a = parameter.preint_sym_expr.subs({"t": tn}).evalf() b = parameter.preint_sym_expr.subs({"t": t}).evalf() if parameter.is_space_dependent: - parameter.dolfin_expression = d.Expression( - sym.printing.ccode((b - a) / dt), degree=3 - ) + c_code = sym.printing.ccode((b - a) / dt) + c_code = c_code.replace("log(", "std::log(") + parameter.dolfin_expression = d.Expression(c_code, degree=3) logger.debug( f"Time-dependent parameter {parameter_name} updated by " f"pre-integrated expression", @@ -2006,7 +2168,9 @@ def dolfin_set_function_values(self, sp, ukey, unew): else: x = d.SpatialCoordinate(sp.compartment.dolfin_mesh) curv = sp.compartment.curv_func - full_expr = d.Expression(sym.printing.ccode(sym_expr), curv=curv, degree=1) + c_code = sym.printing.ccode(sym_expr) + c_code = c_code.replace("log(", "std::log(") + full_expr = d.Expression(c_code, curv=curv, degree=1) ufunc = d.interpolate(full_expr, sp.V) d.assign(sp.u[ukey], ufunc) elif isinstance(unew, d.Expression): diff --git a/smart/model_assembly.py b/smart/model_assembly.py index b8c6ef0..2181ef2 100644 --- a/smart/model_assembly.py +++ b/smart/model_assembly.py @@ -64,6 +64,7 @@ class ParameterType(str, Enum): constant = "constant" expression = "expression" from_xdmf = "from_xdmf" + mesh_quantity = "mesh_quantity" class InvalidObjectException(Exception): @@ -602,6 +603,7 @@ class Parameter(ObjectInstance): is_time_dependent: bool = False is_space_dependent: bool = False compartment: str = "" + mesh_quantity: bool = False def to_dict(self): """Convert to a dict that can be used to recreate the object.""" @@ -731,6 +733,36 @@ def from_xdmf( return parameter + @classmethod + def mesh_quantity( + cls, name, init_val, unit, compartment, group="", notes="", use_preintegration=False + ): + """ " + Initialize as a generic dolfin function over the mesh. + """ + logger.debug(f"Initializing parameter {name} as mesh quantity") + if use_preintegration: + logger.warning( + f"Setting use_preintegration to False for parameter {name}." + "Not currently implemented for parameters given as mesh quantities" + ) + use_preintegration = False + parameter = cls( + name, + init_val, + unit, + group=group, + notes=notes, + use_preintegration=use_preintegration, + ) + parameter.compartment = compartment + # initialize instance + parameter.is_time_dependent = False + parameter.is_space_dependent = True + parameter.type = ParameterType.mesh_quantity + parameter.__post_init__() + return parameter + @classmethod def from_expression( cls, @@ -775,7 +807,9 @@ def from_expression( if is_time_dependent and not is_space_dependent: value = float(sym_expr.subs({"t": 0.0})) else: - dolfin_expression = d.Expression(sym.printing.ccode(sym_expr), t=0.0, degree=3) + c_code = sym.printing.ccode(sym_expr) + c_code = c_code.replace("log(", "std::log(") + dolfin_expression = d.Expression(c_code, t=0.0, degree=3) value = float(sym_expr.subs({"t": 0.0, "x[0]": 0.0, "x[1]": 0.0, "x[2]": 0.0})) parameter = cls( @@ -846,7 +880,7 @@ def __post_init__(self): @property def dolfin_quantity(self): - if self.type == ParameterType.from_xdmf: + if self.type == ParameterType.from_xdmf or self.type == ParameterType.mesh_quantity: return self.dolfin_function * self.unit elif hasattr(self, "dolfin_expression"): return self.dolfin_expression * self.unit @@ -1002,9 +1036,9 @@ def __post_init__(self): f"Creating dolfin object for space-dependent initial condition {self.name}", extra=dict(format_type="log"), ) - self.initial_condition_expression = d.Expression( - sym.printing.ccode(sym_expr), degree=1 - ) + c_code = sym.printing.ccode(sym_expr) + c_code = c_code.replace("log(", "std::log(") + self.initial_condition_expression = d.Expression(c_code, degree=1) elif isinstance(self.initial_condition, Path): pass # keep as path else: @@ -1135,12 +1169,19 @@ class Compartment(ObjectInstance): dimensionality: topological dimensionality (e.g. 3 for volume, 2 for surface) compartment_units: length units for the compartment cell_marker: marker value identifying the compartment in the parent mesh + vel: string expressions for advective velocity field within compartment + deform: string expressions for deformation field within compartment """ name: str dimensionality: int compartment_units: pint.Unit cell_marker: Any + # vel: Union[list[str], list[float]] = [0.0, 0.0, 0.0] + # deform: Union[list[str], list[float]] = [0.0, 0.0, 0.0] + vel: list = dataclasses.field(default_factory=lambda: [0.0, 0.0, 0.0]) + deform: list = dataclasses.field(default_factory=lambda: [0.0, 0.0, 0.0]) + manual_update: bool = False def to_dict(self): "Convert to a dict that can be used to recreate the object." @@ -1170,6 +1211,12 @@ def __post_init__(self): self._usplit = dict() self.V = None self.v = None + self.vel_expr = None + self.vel_func = None + self.deform_expr = None + self.deform_func = None + self.vel_logic = False + self.deform_logic = False def check_validity(self): """ @@ -1743,9 +1790,9 @@ def _post_init_get_integration_measure(self): "volume-surface_to_volume", ]: # intersection of this surface with boundary of source volume(s) - logger.debug( - "DEBUGGING INTEGRATION MEASURE (only fully defined domains are enabled for now)" - ) + # logger.debug( + # "DEBUGGING INTEGRATION MEASURE (only fully defined domains are enabled for now)" + # ) self.measure = self.surface.mesh.dx self.measure_units = self.surface.compartment_units**self.surface.dimensionality @@ -1783,7 +1830,6 @@ def equation_lambda_eval(self, input_type="quantity"): return unit_to_quantity(self._equation_quantity.units) # Seems like setting this as a @property doesn't cause fenics to recompile - @property def form(self): """-1 factor because terms are defined as if they were on the @@ -1802,10 +1848,73 @@ def form(self): mult = u_mask_new else: mult = d.Constant(-1.0, name="-1") + + # alphaExpr = d.Expression("1.0", degree=1) + # Vcur = d.FunctionSpace(self.surface.mesh.dolfin_mesh, "P", 1) + # self.integral_factor = d.interpolate(alphaExpr, Vcur) + # self.integral_factor = alphaExpr + + if self.topology in ["volume", "surface"]: + if self.destination_compartment.deform_logic: + udef = self.destination_compartment.deform_func + Fcur = d.Identity(3) + d.grad(udef) + Jcur = d.det(Fcur) + if self.topology == "surface": + # Nexpr = d.Expression(("x[0]/R", "x[1]/R", "0.0"), degree=1, R=1) + # Vcur = d.VectorFunctionSpace(self.surface.mesh.dolfin_mesh, "P", 1) + # N = d.interpolate(Nexpr, Vcur) + N = self.surface.normals + self.integral_factor = Jcur * d.sqrt( + d.inner(d.dot(N, d.inv(Fcur)), d.dot(N, d.inv(Fcur))) + ) + else: # then volume + self.integral_factor = Jcur + # mult *= Jcur + else: + self.integral_factor = d.Expression("1.0", degree=1) + elif self.topology in [ + "volume_to_surface", + "surface_to_volume", + "volume-volume_to_surface", + "volume-surface_to_volume", + ]: + source_list = list(self.source_compartments.values()) + if ( + self.destination_compartment.deform_logic + and np.all([source.deform_logic for source in source_list]) + and self.surface.deform_logic + ): + # if (self.topology == "volume_to_surface" or + # self.topology == "volume-volume_to_surface"): + # vol_ref = source_list[0] + # else: + # vol_ref = self.destination_compartment + udef = self.surface.deform_func + Fcur = d.Identity(3) + d.grad(udef) + Jcur = d.det(Fcur) + # Nexpr = d.Expression(("x[0]/R", "x[1]/R", "0.0"), degree=1, R=1) + # Vcur = d.VectorFunctionSpace(self.surface.mesh.dolfin_mesh, "P", 1) + # N = d.interpolate(Nexpr, Vcur) + N = self.surface.normals + self.integral_factor = Jcur * d.sqrt( + d.inner(d.dot(N, d.inv(Fcur)), d.dot(N, d.inv(Fcur))) + ) + # mult *= self.integral_factor + elif ( + self.destination_compartment.deform_logic + or np.any([source.deform_logic for source in source_list]) + or self.surface.deform_logic + ): + logger.warning("FIX: Ensure that deformation must be continuous across interface") + self.integral_factor = d.Expression("1.0", degree=1) + else: + self.integral_factor = d.Expression("1.0", degree=1) + if self.axisymm: return ( mult * x[0] + * self.integral_factor * self.equation_lambda_eval(input_type="value") * self.destination_species.v * self.measure @@ -1813,6 +1922,7 @@ def form(self): else: return ( mult + * self.integral_factor * self.equation_lambda_eval(input_type="value") * self.destination_species.v * self.measure diff --git a/smart/solvers.py b/smart/solvers.py index b3b9e24..e488160 100644 --- a/smart/solvers.py +++ b/smart/solvers.py @@ -9,6 +9,8 @@ from .common import Stopwatch from .model_assembly import Compartment +import numpy as np + logger = logging.getLogger(__name__) __all__ = ["smartSNESProblem"] @@ -356,7 +358,14 @@ def assemble_Fnest(self, Fnest): extra=dict(format_type="log"), ) continue - Fvecs[j].append(d.as_backend_type(d.assemble_mixed(form))) + # Fvecs[j].append(d.as_backend_type(d.assemble_mixed(form))) + # assemble_mixed sometimes returns nan (nondeterministic???) + cur_vec = d.assemble_mixed(self.Fforms[j][k]) + count = 1 + while any(np.isnan(cur_vec.get_local())) and count < 10: + cur_vec = d.assemble_mixed(self.Fforms[j][k]) + count += 1 + Fvecs[j].append(d.as_backend_type(cur_vec)) Fj_petsc[j].zeroEntries() for k, vec in enumerate(Fvecs[j]): Fj_petsc[j].axpy(1, vec.vec())