To make it simple to install the package we should consider adding it to conda-forge. However, we should probably make sure that we get a new version of FEniCS released that supports mixed dimensional feature first.

Currently, this feature (see https://github.com/RangamaniLabUCSD/smart/tree/emmetfrancis/multimesh-mpi) only works if every process runs the same number of tasks (for instance, if running 10 conditions, needs to be subdivided into either 1, 2, 5, or 10 processes). Essentially, at the start of each task, the process will hang at dolfin.assemble_mixed within the function, "init_petsc_matnest" inside of smart.solvers until each process reaches that same point, then proceed.

Ideally, we should be able to run any number of processes in parallel as long as we are using the "MPI_SELF" communicator for the mesh and solver, but something within dolfin.assemble_mixed seems to prevent this. As an example that currently fails, consider running "example3_multimeshMPI" (https://github.com/RangamaniLabUCSD/smart/blob/emmetfrancis/multimesh-mpi/examples/example3/example3_multimeshMPI.py) as:

mpiexec -np 3 python3 example3_multimeshMPI.py

The code section of interest reads:

radiusVec = np.logspace(0, 1, num=10)  # currently testing 10 radius values
sim_num = 0
conditions_per_process = int(
    len(radiusVec) / size
)  # currently only works for same number of conditions per process!
ss_vec_cur = np.zeros(conditions_per_process)
for idx in range(rank, len(radiusVec), conditions_per_process):
    curRadius = radiusVec[idx]
    # CODE TO SOLVE THE SYSTEM FOR THIS RADIUS VALUE

In the case of running 3 processes, this results in:

• Process 0 (root): assigned idx = 0, 3, 6, 9
• Process 1: assigned idx = 1, 4, 7
• Process 2: assigned idx = 2, 5, 8

During each run of the loop, there seems to effectively be an MPI barrier at dolfin.assemble_mixed within the function, "init_petsc_matnest" inside of smart.solvers (line 274). This results in process 0 hanging during its fourth run of the loop, and the series of simulations never completes.

Either merge methods into model class or link

Replace main.py with a driver function which takes in a config and set of .jsons

Tables in are not rendered as expected when logged (see e.g https://rangamanilabucsd.github.io/smart/examples/example3/example3.html and screenshot).

One way to fix this is to add a newline after logging the level. However, we probably don't want to make this the default behaviour, so it would perhaps make sense to add a new FormatType for table (e.g in

The MeshView.cell_map() and vertex_map() methods have a value of 0 when there is no connection - however 0 could be the index of a cell or vertex that is connected - so we need to check for this. Also, there is a message that pops up when a connection can't be found (which is very often when a surface mesh is a superset of the boundary of a volume mesh), this leads to thousands of messages which fill the screen and make it difficult to see what stubs is doing.

BoundingBoxTree() MPI bug (#38)

These will be fixed in the forked repository:
https://github.com/justinlaughlin/dolfinstubs/issues/new

According to the joss review checklist we need to add information about how to contribute. A common way to do this is to add a CONTRIBUTING.md, see e.g this file in simcardems

Every time we make a release, we would like to also upload the new version to pypi. This can be done very easy with the following action, see e.g https://github.com/ComputationalPhysiology/simcardems/blob/main/.github/workflows/main.yml

Currently models can only be setup in 3d with 2d boundary conditions. The ability to setup the system in 0d (i.e. well-mixed) and solve via e.g. sympy solve_ivp would allow quicker testing of model parameters

Currently there is no support for specifying Dirichlet boundary conditions.

Some things when left out of the config file, e.g. plot styles, can lead to exceptions. If we set default values to be overwritten this can be prevented.

Other values which are conditionally present should be explicitly checked if it's set.

When a subsurface is a superset of the boundary of a subvolume, we need to use the "nan" (max size_t) values in cell_map to restrict the mesh.

Consider a volume V:

• function u,
• boundary measure ds
  A surface adjacent to the volume S:
• function u_
• measure dx_

We can't assemble u*u_*ds because UFL raises the exception "codim != 0 or 1 - Not (yet) implemented". Even if we restrict ds(z) where z is the marker value of S, this will not work.

We can assemble u*u_*dx_, but this results in an erroneous result if (S intersect V) != S

Need to map "volume facet meshfunction" values to "surface cell meshfunction" values so that we can select a specific portion of the measure dx_. Changes are in stubs/mesh.py

Currently unit conversions assume that each surface only has one associated volume. 2d species need to have a defined "parent compartment" so the correct scaling factor is applied.

https://fenicsproject.org/qa/2532/how-to-set-up-interior-neumann-boundaries/

All of these have some sort of tabular structure - exporting the data into a LaTeX format would be an excellent feature for user-experience and publication.

Currently equation_lambda_eval() will properly satisfy total unit conversions (ie conversion of a flux's units to the type required by the PDE), but it does not ensure that units are converted on the inside of a flux (e.g. u+k where u and k are consistent but different units). The solution is to simply input the pint quantity into equation_lambda but certain expressions don't work (I believe some expressions with exponents)

model._sorted_compartments is naively computed (num_dofs), but more importantly it can lead to inconsistent compartment ordering across processors. Probably best to keep it in serial, and modify it to just exclude non-active compartments in parallel.

Dear authors,

The software installation routine is straightforward, and there is no special issue with it. However, I found the package a bit difficult to get started with, especially for a normal user. If the package is supposed to be used by advanced users only (which I think is not), it should be explicitly mentioned.

• After the installation using the docker approach, how is the user supposed to work with the examples? I think there is no way in the current setup. Jupyter is not installed in the docker image, and there is no port forwarding on the docker run command. I would say a non-advanced user can never find out how to do that, so it can be beneficial to guide the users on how to run the example codes and see how the library works in action. I did it by running jupyter notebook --ip 0.0.0.0 --port 8888 --no-browser --allow-root inside the container (after forwarding the port in the docker execution command), in which the inline PyVista visualization didn't still work (which is because PyVista is not installed inside the container).
• As a follow-up to the above point, it can be nice if the user can have access to the example codes in separate files without the Jupyter notebooks, similar to example 3. These codes (accompanied by some documentation on how to run them) can show a normal user how to use the package after the installation.

P.S. This issue is related to openjournals/joss-reviews#5580

Probably will require a fork of dolfin.

Detailed explanation of issue here:
https://fenicsproject.discourse.group/t/assembly-on-mixed-domain-hangs-segfaults-in-parallel-dependent-on-partitioning/7583/2

Before we can submit the package to joss we need to write the software paper. I think the material in the overleaf document is more than enough, but we need to convert it to the correct format.

Alternative names for stubs has been suggested. One alternative is

Spatial Modelling Algorithms for Reaction-Transport [systems|models|equations]

with the acronym SMART. We need to

• agree on an updated name and acronym
• update the name in the file tree and documentation

Would be good to do this sooner rather than later. @emmetfrancis Are you all settled on the name SMART?

Currently stubs automatically generates plots for min, mean, and max taken over the entire domain. a fenics function, u, can be evaluated at a point with the syntax u(x,y,z). The functionality to generate plots over time for specific coordinates should be implemented.

Looks like we can't solve mixed-dimension problems using dolfin's built in solvers:

dolfin.fem.problem.MixedNonlinearVariationalProblem() (Jacobian must be provided)
d.fem.solving.solve(F==0,u) (which then calls _solve_varproblem(), then MixedNonlinearVariationalProblem(), however this leads to an assertion error of an incorrectly sized Jacobian [generated in lines 286-291 of dolfin.fem.solving])
dolfin.fem.solving.assemble_mixed_system() (raises RuntimeError "assemble_mixed_system is not available for non-linear problems yet")

After some investigation, it looks like MixedNonlinearVariationalProblem() doesn't construct the Jacobian correctly for vector function spaces in a mixed function spaces.

Since MixedLinearVariationalProblem() seems to be working, the most straightforward solution seems to be the following:

1. Manually construct the Jacobian for our problem (in particular vector function spaces)
2. See if we can use this Jacobian in Newton's method with MixedNonlinearVariationalProblem/Solver (SNES solver does not work yet with these methods)
3. If 2) does not work, then formulate the problem as a linearized one and manually wrap our own Newton's method.

Dear authors,

The software has standard and comprehensive automated tests implemented using GitHub workflows. However, it's not very clear how a user/developer can test a local installation of the code. Inside the docker container, running pytest (or python3 -m pytest) command returns this error:

_main__.py: error: unrecognized arguments: --cov=smart --cov-report --cov-report term-missing
  inifile: /home/shared/pytest.ini
  rootdir: /home/shared

P.S. This issue is related to openjournals/joss-reviews#5580

Should be a function of:

• "multiphysics" solver (iterative, monolithic, ...)
• nonlinear solver (newton, picard)
• linear solver (most likely with krylov settings to pass to dolfin)

Purpose is to separate from the Model class to improve modularity

One challenge with some om this classes is the attributes are set dynamically during initialisation. For example, when we create a Species is noes not know about the compartment, but it has a method called vscalar which look as follows

class Species(ObjectInstance):
    ....
    @cached_property
    def vscalar(self):
        return d.TestFunction(sub(self.compartment.V, 0, True))

This is a problem because an instance of Species don't have an attribute compartment, but this is set dynamically during the model initialisation, see https://github.com/RangamaniLabUCSD/stubs/blob/development/stubs/model.py#L582

One way to make this a bit more explicit it is initialise compartment to a default values (e.g None), and raise some error indicating that the object is not complete, e.g

@dataclass
class Species(ObjectInstance):
    compartment: Optional[Compartment] = None
    ....
    @cached_property
    def vscalar(self):
        if self.compartment is None:
            raise RunTimeError("Species is not complete")
        return d.TestFunction(sub(self.compartment.V, 0, True))

Ideally, we should rethink how we construct these object and try to avoid to tight coupling between them. Does Species really need to know the compartments, or is it enough that the model has this information?

We could probably get some inspiration from the Builder pattern https://refactoring.guru/design-patterns/builder which tries to address a similar issue.

Allow user to define expressions in config/model file. Seems like the cleanest route for implementing boolean logic/unit-step functions e.g.
f = 5 + 2*Heaviside(u-1)
can be written as
f=Expression('u < 1 ? 5 : 7',u=u)

For some reason, the logs for warnings starts with .... (see https://rangamanilabucsd.github.io/smart/examples/example3/example3.html and screenshot)

I guess this is related to the filler_char in https://github.com/RangamaniLabUCSD/smart/blob/97ad37992f8cf07d5ee53fbb1491adfc0712b6e6/smart/config.py#LL132C18-L132C18

Is this something we want?

Eventually, we want to explore the possibility of converting our models that have time-dependent parameters to SBML. This doesn't seem to be supported yet, and we have to be able to convert both functions and CSV files to SBML.

One possibility we can explore is using Rules (see the spec here) to express our equations and define time as a standalone parameter, but this should be saved for after we have our base SBML capabilities in place.

Of course, native support is always preferred, and Chris suggested raising an issue to the SBML team to see if they can add support if we can instead of doing some workaround.

In several of the examples (e.g https://rangamanilabucsd.github.io/smart/examples/example3/example3.html) we see a lot of warnings from pandas

usr/lib/python3/dist-packages/pandas/core/dtypes/cast.py:1579: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  result[:] = values

We should try to figure out why we get these warnings and fix it.

Pre-integration allows efficient computation when purely time-dependent functions are very stiff, however sympy takes a very long time to integrate certain functions e.g. Heaviside. For now, users can provide a pre-defined indefinite integral of the time-dependent function and use the is_preintegrated flag.

Current issues with using MPI to parallelize a mesh in SMART:

• Solver will not work unless there is at least one instance of each compartment within each partition of the mesh.
• A partition must not have any ghost faces, i.e. any cell-face-cell interface must have both cells owned by the same process

An example illustrating these problems will be added soon.

Compute flux balance at boundaries. User should define which parameters are free. Export system of ODEs and solve for roots.

Currently, this project is using readthedocs as a backend for its documentation.
If we want a tighter integration with Github, we could use JupyterBook, as for instance done for
https://github.com/ComputationalPhysiology/simcardems
https://computationalphysiology.github.io/simcardems/
and
https://github.com/jorgensd/dolfinx-tutorial
https://jsdokken.com/dolfinx-tutorial/

It makes it easier to test locally. @emmetfrancis @justinlaughlin @finsberg any thoughts?

It could possibly be a useful feature to report parameters that are unused by any model. At the least, it could reduce the database size. This could be useful if folks eventually develop databases of reactions and try to combine/merge them.

Some sort of wrapper class or function to perform parameter sweeps and analyze the results.

Common values that a user may want to vary in a parameter sweep:

• parameters
• initial concentrations
• diffusion coefficients
• configuration/solver settings (e.g. different meshes, solver tolerances, etc.)

A wrapper function to load the data from multiple simulations and compare them would be a reasonable next step. Data from simulations (not the raw 3d time-series data but rather extracted values such as min, mean, max, etc) is serialized using pickle. Data from multiple pickles can be loaded and plotted against each other.

Long term consideration...possible integration with PSUADE? https://computing.llnl.gov/projects/psuade/software

The code is currently unused
https://github.com/search?q=repo%3ARangamaniLabUCSD%2Fsmart%20MixedAssemblerTemp&type=code

@emmetfrancis do you know if it is still needed?
The reason (according to the comments) was:

Our custom assembler (something about dolfin's init_global_tensor was not

    # correct so we manually initialize the petsc matrix and then wrap with dolfin)
    # This assembly routine is the exact same as d.assemble_mixed() except
    # init_global_tensor() is commented out
    # Thanks to Prof. Kamensky and his student for the idea
    # https://github.com/hanzhao2020/PENGoLINS/blob/main/PENGoLINS/cpp/transfer_matrix.cpp

FEniCS doesn't support parallelization on SubMeshes. The current workaround solution simply gives the entire submesh to each CPU but this is grossly inefficient when a submesh has >1e5 dofs. Two options are:

1. Modify cbcpost code (python 2) and implement.
2. Project submesh onto its parent mesh (values are 0 outside of the submesh - potentially viable if submesh contains 50+% of the vertices of the parent)

Modelers may be interested in tracking a scalar variable which depends on e.g. reactive fluxes. Example: Ca2+ concentration in a single ER store assuming fast diffusion.

Implement ability to generate a reaction network from a model file.

Pysb and NetworkX are two packages to check out.

Currently stubs uses the fancy_print function to log everything. This is essentially a wrapper around the print function for creating coloured output as well as some options for adding new lines and titles. You can also provide a filename in order to dump the output to a file.

The drawback with using this is that there is now way to customise the amount of data that is logged. This is usually done using a logger which provides a way to log information at different log levels (see e.g https://realpython.com/python-logging/ for a gentle introduction).

I have now been looking into how we can swap out the fancy_print function with proper logging using the logging module, but before implementing it I would like to get some feedback on it.

The loggers in the logging module consists of handlers (e.g whether you want output to a file or to the console), and each handler has a formatter where you can customise the way the output is logged.

Consider the following example where I try to use a custom formatter to mimic some of the features from the fancy_print function

import logging
import dolfin
from termcolor import colored
from stubs.common import _fancy_print as fancy_print

comm = dolfin.MPI.comm_world
rank = comm.rank
size = comm.size

base_format = (
    "%(asctime)s %(rank)s%(name)s - %(levelname)s %(message)s (%(filename)s:%(lineno)d)"
)

class StubsFormatter(logging.Formatter):
    def format(self, record: logging.LogRecord) -> str:

        record.rank = f"CPU {rank}: " if size > 1 else ""
        format_type = getattr(record, "format_type", "default")
        if format_type == "title":
            record.msg = f" {record.msg} ".center(40, "=")
        color = getattr(record, "color", "white")
        log_fmt = colored(base_format, color)
        formatter = logging.Formatter(log_fmt)

        out = formatter.format(record)

        new_lines = getattr(record, "new_lines", (0, 0))
        prefix = "\n" * new_lines[0]
        postfix = "\n" * new_lines[1]
        return prefix + out + postfix



stream_handler = logging.StreamHandler()
stream_handler.setFormatter(StubsFormatter())

file_handler = logging.FileHandler("file.log")
file_handler.setFormatter(StubsFormatter())

logger = logging.getLogger(__name__)
logger.addHandler(stream_handler)
logger.addHandler(file_handler)

logger.setLevel(logging.DEBUG)


logger.info("Base text with logger")
fancy_print("Base text with fancy_print")

logger.info("Logger with color", extra=dict(color="magenta"))
fancy_print("Fancy with color", text_color="magenta")

logger.info("Title text with logger", extra=dict(format_type="title"))
fancy_print("Title text with fancy_print", format_type="title")

logger.info("Base text with logger with newline before", extra=dict(new_lines=(1, 0)))
fancy_print("Base text with fancy_print with newline before", new_lines=(1, 0))

Here I am creating two different handlers, one that outputs the logs to a file called file.log and one that output the same logs to the console. Both of these handlers gets a custom formatter which will mimic the fancy_print function. The custom formatter class (StubsFormatter) has a method called format which are called automatically when you log a message and which operates on an object of type LogRecord. You can pass extra arguments to this LogRecord object using the extra keyword argument which is one way that we could pass on information about how the message should be logged.

Currently I have implemented a way to add a color (using termcolor which is currently used in fancy_print), adding new lines before and after the message and creating a title by adding = on both sides of the text. I also added the CPU rank in case we are running on multiple processes.

Would be great to get some feedback on whether this is a good solution or not, or if you have any other suggestions?

Translate our models into SBML and vice-versa

Dear authors,

Generally speaking, the software runs as described in the documentation. However, there are very minor comments in this regard:

• The software claims to be high-performance, but no performance test is reported. Having a small test on example 3 (as the relevant example on parallelization) can be nice in this regard.
• Is there any way to reproduce Figure 1 of the paper? Example 2 performs the simulation on the mentioned model on a simple geometry.

P.S. This issue is related to openjournals/joss-reviews#5580

Related to #14

Similar to here we should create and publish a new docker image for every release. We should also make sure to support both arm and amd architecture, but it should be more or less copy-paste from here

mesh.mf_map[(a,b)] will contain a mesh function with the intersection of mesh with mesh (id=a) and mesh (id=b)

and dx_map will be the Measure() used to integrate

Saving solutions out as hdf5 allows them to be loaded back in as dolfin functions. This is important for scenarios such as:

• Solving for steady-state (which could be spatially heterogeneous) and then using that steady-state as initial conditions
• Recovering from a failed run with smaller time-steps

Dear authors,

The documentation of the software is nicely done and contains an overview of the underlying mathematical and computational aspects too. However, you may consider improving it with these suggestions:

• It can be beneficial to talk a little bit about the visualization of results in an external postprocessor like ParaView. For all the examples, the results are written on the disc, but the documentation and examples don’t mention how to deal with them except by using the provided analysis tools. If you assume that the user should already know about how FEniCS and its IO pipelines work, you need to mention these prerequisites explicitly (also in the paper).
• For the examples, it can be beneficial to mention where the values of the model parameters are coming from. This can help the users build their own models from scratch. For example, are the chosen values of parameters in example 6 based on the model of Doi et al.? If not, how are the values chosen?
• It can be beneficial to mention that some of the cells of the Jupyter files take longer to compute, maybe with an approximation of the time that the user needs to wait for the results to come. For instance, the loop in example 3 seems to never end (it took half an hour to finish for me I think). The same problem exists with example 6.
• How is an advanced user supposed to change the performance-related configs, like the default PETSc preconditioner and solver types and parameters used by SMART? Is that a possibility via the configCur.solver.update function?

P.S. This issue is related to openjournals/joss-reviews#5580