From 252e32f97edd067b0284d8ffead08ac61cabc08f Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Tue, 9 Apr 2024 23:03:14 -0700 Subject: [PATCH 1/7] Convert equation strings to equations, fix nan returned sometimes during solving --- examples/example3/example3_withaxisymm.ipynb | 4 +- smart/model.py | 6 +++ smart/model_assembly.py | 54 ++++++++++---------- smart/solvers.py | 7 ++- 4 files changed, 40 insertions(+), 31 deletions(-) diff --git a/examples/example3/example3_withaxisymm.ipynb b/examples/example3/example3_withaxisymm.ipynb index 93b402db..f57e54a0 100644 --- a/examples/example3/example3_withaxisymm.ipynb +++ b/examples/example3/example3_withaxisymm.ipynb @@ -367,7 +367,7 @@ "outputs": [], "source": [ "# Base mesh\n", - "domain, facet_markers, cell_markers = mesh_tools.create_2Dcell(outerExpr=f\"r**2 + (z-({curRadius}+1))**2 - {curRadius}**2\", innerExpr=\"\", hEdge=0.1)\n", + "domain, facet_markers, cell_markers = mesh_tools.create_2Dcell(outerExpr=f\"r**2 + (z-({curRadius}+1))**2 - {curRadius}**2\", innerExpr=\"\", hEdge=0.025, hInnerEdge=0.025)\n", "visualization.plot_dolfin_mesh(domain, cell_markers)" ] }, @@ -425,7 +425,7 @@ " \"final_t\": 1,\n", " \"initial_dt\": 0.01,\n", " \"time_precision\": 6,\n", - " \"attempt_timestep_restart_on_divergence\": True\n", + " # \"attempt_timestep_restart_on_divergence\": True\n", " }\n", ")\n", "config_cur.flags.update({\"axisymmetric_model\": True})" diff --git a/smart/model.py b/smart/model.py index dc532f19..7f6baab5 100644 --- a/smart/model.py +++ b/smart/model.py @@ -472,6 +472,12 @@ def _init_2_1_reactions_to_symbolic_strings(self): raise ValueError( "Reaction %s does not seem to have an associated equation" % reaction.name ) + if reaction.eqn_f_str == "": + reaction.eqn_str = f"-{reaction.eqn_r_str}" + elif reaction.eqn_r_str == "": + reaction.eqn_str = f"{reaction.eqn_f_str}" + else: + reaction.eqn_str = f"{reaction.eqn_f_str}-{reaction.eqn_r_str}" def _init_2_2_check_reaction_validity(self): """Confirms that all reactions have parameters/species defined""" diff --git a/smart/model_assembly.py b/smart/model_assembly.py index d7a4651d..ecf0bd27 100644 --- a/smart/model_assembly.py +++ b/smart/model_assembly.py @@ -242,7 +242,29 @@ def print_to_latex( for row in range(df.shape[0]): for col in range(df.shape[1]): - if isinstance(df.iat[row, col], str): + if df.columns[col] == "eqn_str": + cur_str = df.iat[row, col] + cur_str = "$" + cur_str + "$" + idx = 0 + while idx < len(cur_str): + if cur_str[idx] == "_": + cur_str = cur_str[0 : idx + 1] + "{" + cur_str[idx + 1 :] + isMath = False + testIdx = idx + 2 + while not isMath: + if not cur_str[testIdx].isalnum(): + if cur_str[testIdx] == "*": + cur_str = cur_str[0:testIdx] + "} " + cur_str[testIdx + 1 :] + else: + cur_str = cur_str[0:testIdx] + "}" + cur_str[testIdx:] + isMath = True + else: + testIdx += 1 + idx = testIdx + 2 + else: + idx += 1 + df.iloc[row, col] = cur_str + elif isinstance(df.iat[row, col], str): cur_str = df.iat[row, col] if "_" in cur_str: df.iloc[row, col] = cur_str.replace("_", "\_") @@ -776,31 +798,6 @@ def print( s.latex_name super().print(tablefmt, self.properties_to_print, filename, max_col_width) - def print_to_latex( - self, - properties_to_print=None, - max_col_width=None, - sig_figs=2, - return_df=False, - ): - # properties_to_print = ["_latex_name"] - # properties_to_print.extend(self.properties_to_print) - df = super().print_to_latex(properties_to_print, max_col_width, sig_figs, return_df=True) - # fix dof_index - for col in df.columns: - if col == "dof_index": - df[col] = df[col].astype(int) - # fix name - # get the column of df that contains the name - # this can be more robust - # df.columns - - if return_df: - return df - else: - with pandas.option_context("max_colwidth", 1000): - logger.info(df.to_latex(escape=False, longtable=True, index=False)) - @dataclass class Species(ObjectInstance): @@ -841,7 +838,6 @@ class Species(ObjectInstance): def to_dict(self): "Convert to a dict that can be used to recreate the object." keys_to_keep = [ - "name", "initial_condition", "concentration_units", "D", @@ -1139,7 +1135,7 @@ class ReactionContainer(ObjectContainer): def __init__(self): super().__init__(Reaction) - self.properties_to_print = ["lhs", "rhs", "eqn_f_str", "eqn_r_str"] + self.properties_to_print = ["lhs", "rhs", "eqn_str"] def print( self, @@ -1217,6 +1213,7 @@ class Reaction(ObjectInstance): track_value: bool = False eqn_f_str: str = "" eqn_r_str: str = "" + eqn_str: str = "" group: str = "" axisymm: bool = False @@ -1234,6 +1231,7 @@ def to_dict(self): "track_value", "eqn_f_str", "eqn_r_str", + "eqn_str", "group", "axisymm", ] diff --git a/smart/solvers.py b/smart/solvers.py index fe124ba4..a9f6be2e 100644 --- a/smart/solvers.py +++ b/smart/solvers.py @@ -8,6 +8,7 @@ from .common import Stopwatch from .model_assembly import Compartment +import numpy as np logger = logging.getLogger(__name__) @@ -350,7 +351,11 @@ def assemble_Fnest(self, Fnest): Fvecs.append([]) for k in range(len(self.Fforms[j])): # , tensor=d.PETScVector(Fvecs[idx])) - Fvecs[j].append(d.as_backend_type(d.assemble_mixed(self.Fforms[j][k]))) + cur_vec = d.assemble_mixed(self.Fforms[j][k]) + # assemble_mixed sometimes returns nan (nondeterministic???) + while any(np.isnan(cur_vec.get_local())): + cur_vec = d.assemble_mixed(self.Fforms[j][k]) + Fvecs[j].append(d.as_backend_type(cur_vec)) # TODO: could probably speed this up by not using axpy if there is only one subform # sum the vectors Fj_petsc[j].zeroEntries() From d0cffe8edc269122d8629a54651b4b4b70922ba8 Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Wed, 17 Apr 2024 10:58:49 -0700 Subject: [PATCH 2/7] Fixes to formatting math in reactions tables --- smart/model_assembly.py | 5 +++- smart/test_nan.ipynb | 57 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 smart/test_nan.ipynb diff --git a/smart/model_assembly.py b/smart/model_assembly.py index ecf0bd27..d41316af 100644 --- a/smart/model_assembly.py +++ b/smart/model_assembly.py @@ -252,7 +252,7 @@ def print_to_latex( isMath = False testIdx = idx + 2 while not isMath: - if not cur_str[testIdx].isalnum(): + if not (cur_str[testIdx].isalnum() or cur_str[testIdx] == "_"): if cur_str[testIdx] == "*": cur_str = cur_str[0:testIdx] + "} " + cur_str[testIdx + 1 :] else: @@ -261,6 +261,9 @@ def print_to_latex( else: testIdx += 1 idx = testIdx + 2 + elif cur_str[idx] == "*": + cur_str = cur_str[0:idx] + " " + cur_str[idx + 1 :] + idx += 1 else: idx += 1 df.iloc[row, col] = cur_str diff --git a/smart/test_nan.ipynb b/smart/test_nan.ipynb new file mode 100644 index 00000000..3cdf236c --- /dev/null +++ b/smart/test_nan.ipynb @@ -0,0 +1,57 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from smart import mesh_tools\n", + "import dolfin as d\n", + "import numpy as np\n", + "curRadius = 1.0\n", + "domain, facet_markers, cell_markers = mesh_tools.create_2Dcell(outerExpr=f\"r**2 + (z-({curRadius}+1))**2 - {curRadius}**2\", innerExpr=\"\", hEdge=0.01, hInnerEdge=0.01)\n", + "cyto_mesh = d.create_meshview(cell_markers, 1)\n", + "pm_mesh = d.create_meshview(facet_markers, 10)\n", + "cyto_func_space = d.FunctionSpace(cyto_mesh, \"P\", 1)\n", + "test_expr = d.Expression(f\"1/(pow(x[0],2)+pow(x[1],2)+pow(x[2]-({curRadius}+1),2)+0.1)\", degree=1)\n", + "cyto_func = d.interpolate(test_expr, cyto_func_space)\n", + "xCyto = d.SpatialCoordinate(cyto_mesh)\n", + "dxCyto = d.Measure(\"dx\", cyto_mesh)\n", + "testFunc = d.TestFunction(cyto_func_space)\n", + "Nmax = 1e4\n", + "N = 0\n", + "while N < Nmax:\n", + " cur_val = d.assemble_mixed(xCyto[0]*cyto_func*testFunc*dxCyto)\n", + " if np.any(np.isnan(cur_val.get_local())):\n", + " print(f\"Current integral is {cur_val.get_local()}\")\n", + " break\n", + " N += 1\n", + " print(N)\n", + "\n", + "print(f\"Current integral is {cur_val.get_local()}\")" + ] + } + ], + "metadata": { + "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": 2 +} From 23e5b5d3481728dd32f0fb40bac71681074968fb Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Wed, 17 Apr 2024 11:00:23 -0700 Subject: [PATCH 3/7] Move nan test --- smart/test_nan.ipynb | 57 -------------------------------------------- 1 file changed, 57 deletions(-) delete mode 100644 smart/test_nan.ipynb diff --git a/smart/test_nan.ipynb b/smart/test_nan.ipynb deleted file mode 100644 index 3cdf236c..00000000 --- a/smart/test_nan.ipynb +++ /dev/null @@ -1,57 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from smart import mesh_tools\n", - "import dolfin as d\n", - "import numpy as np\n", - "curRadius = 1.0\n", - "domain, facet_markers, cell_markers = mesh_tools.create_2Dcell(outerExpr=f\"r**2 + (z-({curRadius}+1))**2 - {curRadius}**2\", innerExpr=\"\", hEdge=0.01, hInnerEdge=0.01)\n", - "cyto_mesh = d.create_meshview(cell_markers, 1)\n", - "pm_mesh = d.create_meshview(facet_markers, 10)\n", - "cyto_func_space = d.FunctionSpace(cyto_mesh, \"P\", 1)\n", - "test_expr = d.Expression(f\"1/(pow(x[0],2)+pow(x[1],2)+pow(x[2]-({curRadius}+1),2)+0.1)\", degree=1)\n", - "cyto_func = d.interpolate(test_expr, cyto_func_space)\n", - "xCyto = d.SpatialCoordinate(cyto_mesh)\n", - "dxCyto = d.Measure(\"dx\", cyto_mesh)\n", - "testFunc = d.TestFunction(cyto_func_space)\n", - "Nmax = 1e4\n", - "N = 0\n", - "while N < Nmax:\n", - " cur_val = d.assemble_mixed(xCyto[0]*cyto_func*testFunc*dxCyto)\n", - " if np.any(np.isnan(cur_val.get_local())):\n", - " print(f\"Current integral is {cur_val.get_local()}\")\n", - " break\n", - " N += 1\n", - " print(N)\n", - "\n", - "print(f\"Current integral is {cur_val.get_local()}\")" - ] - } - ], - "metadata": { - "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": 2 -} From d11146443f3ca4e2c00c1aa5b291616a21b4c8ad Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Tue, 30 Apr 2024 08:19:10 -0700 Subject: [PATCH 4/7] Fixes to formatting latex tables for pint quantities --- smart/model_assembly.py | 45 +++++++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/smart/model_assembly.py b/smart/model_assembly.py index 1ac1c38b..ae163154 100644 --- a/smart/model_assembly.py +++ b/smart/model_assembly.py @@ -219,8 +219,8 @@ def print_to_latex( self, properties_to_print=None, max_col_width=None, - sig_figs=2, - return_df=True, + sig_figs=3, + return_df=False, ): """ Print object properties in latex format. @@ -245,6 +245,7 @@ def print_to_latex( if df.columns[col] == "eqn_str": cur_str = df.iat[row, col] cur_str = "$" + cur_str + "$" + cur_str = cur_str.replace("**", "^") idx = 0 while idx < len(cur_str): if cur_str[idx] == "_": @@ -284,15 +285,15 @@ def print_to_latex( new_list[i] = "`" cur_str = "".join(new_list) df.iloc[row, col] = cur_str - for col in df.columns: - # Convert quantity objects to unit - if isinstance(df[col].iloc[0], pint.Quantity): - # if tablefmt=='latex': - df[col] = df[col].apply(lambda x: f"${x:0.{sig_figs}e~Lx}$") - - if col == "idx": - df = df.drop("idx", axis=1) + # Convert quantity objects to unit + elif isinstance(df.iat[row, col], pint.Quantity): + x = df.iat[row, col] + if isinstance(x.magnitude, str): + df.iloc[row, col] = f"{x:s~P}" + else: + df.iloc[row, col] = f"${x:0.{sig_figs}e~Lx}$" + for col in df.columns: if "_" in col: df = df.rename(columns={col: col.replace("_", "\_")}) @@ -300,7 +301,7 @@ def print_to_latex( return df else: with pandas.option_context("max_colwidth", 1000): - logger.info(df.to_latex(escape=False, longtable=True, index=False)) + logger.info(df.to_latex(escape=False, longtable=True, index=True)) def get_pandas_dataframe_formatted( self, @@ -351,11 +352,14 @@ def print( ) # # Change certain df entries to best format for printing - for col in df.columns: - # Convert quantity objects to unit - if isinstance(df[col].iloc[0], pint.Quantity): - # if tablefmt=='latex': - df[col] = df[col].apply(lambda x: f"{x:0.{sig_figs}e~P}") + for row in range(df.shape[0]): + for col in range(df.shape[1]): + if isinstance(df.iat[row, col], pint.Quantity): + x = df.iat[row, col] + if isinstance(x.magnitude, str): + df.iloc[row, col] = f"{x:s~P}" + else: + df.iloc[row, col] = f"{x:0.{sig_figs}e~P}" # print to file if filename is None: @@ -472,7 +476,7 @@ def __init__(self): self.properties_to_print = [ "_quantity", - "is_time_dependent", + # "is_time_dependent", "sym_expr", "notes", "group", @@ -785,8 +789,7 @@ def __init__(self): self.properties_to_print = [ "compartment_name", "_Diffusion", - "initial_condition", - "concentration_units", + "_Initial_Concentration", ] def print( @@ -799,6 +802,7 @@ def print( for s in self: s.D_quantity s.latex_name + s.initial_condition_quantity super().print(tablefmt, self.properties_to_print, filename, max_col_width) @@ -1138,7 +1142,7 @@ class ReactionContainer(ObjectContainer): def __init__(self): super().__init__(Reaction) - self.properties_to_print = ["lhs", "rhs", "eqn_str"] + self.properties_to_print = ["lhs", "rhs", "eqn_str", "topology"] def print( self, @@ -1217,6 +1221,7 @@ class Reaction(ObjectInstance): eqn_f_str: str = "" eqn_r_str: str = "" eqn_str: str = "" + topology: str = "" group: str = "" axisymm: bool = False has_subdomain: bool = False From 3deb8644a3ea9b2f43b12b4e743895dd28022f68 Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Tue, 30 Apr 2024 08:24:46 -0700 Subject: [PATCH 5/7] Revert changes to example3withaxisymm --- examples/example3/example3_withaxisymm.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/example3/example3_withaxisymm.ipynb b/examples/example3/example3_withaxisymm.ipynb index f57e54a0..d91a4fca 100644 --- a/examples/example3/example3_withaxisymm.ipynb +++ b/examples/example3/example3_withaxisymm.ipynb @@ -367,7 +367,7 @@ "outputs": [], "source": [ "# Base mesh\n", - "domain, facet_markers, cell_markers = mesh_tools.create_2Dcell(outerExpr=f\"r**2 + (z-({curRadius}+1))**2 - {curRadius}**2\", innerExpr=\"\", hEdge=0.025, hInnerEdge=0.025)\n", + "domain, facet_markers, cell_markers = mesh_tools.create_2Dcell(outerExpr=f\"r**2 + (z-({curRadius}+1))**2 - {curRadius}**2\", innerExpr=\"\", hEdge=0.1)\n", "visualization.plot_dolfin_mesh(domain, cell_markers)" ] }, @@ -627,7 +627,7 @@ " \"final_t\": 1,\n", " \"initial_dt\": 0.01,\n", " \"time_precision\": 6,\n", - " \"attempt_timestep_restart_on_divergence\": True,\n", + " # \"attempt_timestep_restart_on_divergence\": True,\n", " }\n", " )\n", " config_cur.flags.update({\"axisymmetric_model\": True})\n", From c35a3a1826ad83b77472f9539156f2b05b52a95c Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Tue, 7 May 2024 22:35:15 -0700 Subject: [PATCH 6/7] Change back setting vector values in restrict function for species --- smart/mesh_tools.py | 9 ++++++--- smart/model.py | 6 +++++- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py index 6d4c2ec7..576a1be7 100644 --- a/smart/mesh_tools.py +++ b/smart/mesh_tools.py @@ -1438,16 +1438,19 @@ def load_mesh( if extra_keys is None: extra_keys = [] - if not pathlib.Path(filename).is_file(): + if isinstance(filename, str): + filename = pathlib.Path(filename) + + if not filename.is_file(): raise FileNotFoundError(f"File {filename} does not exists") if mesh is None: mesh = d.Mesh(comm) - with d.HDF5File(comm, pathlib.Path(filename).with_suffix(".h5").as_posix(), "r") as hdf5: + with d.HDF5File(comm, filename.with_suffix(".h5").as_posix(), "r") as hdf5: hdf5.read(mesh, "/mesh", False) - dim = mesh.geometric_dimension() + dim = mesh.topology().dim() # geometric_dimension() mf_cell = d.MeshFunction("size_t", mesh, dim) mf_facet = d.MeshFunction("size_t", mesh, dim - 1) subdomains = [] diff --git a/smart/model.py b/smart/model.py index f6993b25..e11877bc 100644 --- a/smart/model.py +++ b/smart/model.py @@ -992,7 +992,11 @@ def _init_4_7_set_initial_conditions(self): # restrict to specified subdomain u_cur = self.cc[species.compartment_name].u[ukey] u_new = create_restriction(u_cur, species.subdomain_data, species.subdomain_val) - u_cur.assign(u_new) + values = u_cur.vector().get_local() + values_new = u_new.vector().get_local() + values[species.dof_map] = values_new[species.dof_map] + u_cur.vector().set_local(values) + u_cur.vector().apply("insert") def _init_5_1_reactions_to_fluxes(self): """Convert reactions to flux objects""" From 75bbd762eade527ba9abc5f747cbfcaf1220f0ee Mon Sep 17 00:00:00 2001 From: emmetfrancis Date: Tue, 7 May 2024 23:16:56 -0700 Subject: [PATCH 7/7] Remove nan workaround --- smart/solvers.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/smart/solvers.py b/smart/solvers.py index a9f6be2e..fe124ba4 100644 --- a/smart/solvers.py +++ b/smart/solvers.py @@ -8,7 +8,6 @@ from .common import Stopwatch from .model_assembly import Compartment -import numpy as np logger = logging.getLogger(__name__) @@ -351,11 +350,7 @@ def assemble_Fnest(self, Fnest): Fvecs.append([]) for k in range(len(self.Fforms[j])): # , tensor=d.PETScVector(Fvecs[idx])) - cur_vec = d.assemble_mixed(self.Fforms[j][k]) - # assemble_mixed sometimes returns nan (nondeterministic???) - while any(np.isnan(cur_vec.get_local())): - cur_vec = d.assemble_mixed(self.Fforms[j][k]) - Fvecs[j].append(d.as_backend_type(cur_vec)) + Fvecs[j].append(d.as_backend_type(d.assemble_mixed(self.Fforms[j][k]))) # TODO: could probably speed this up by not using axpy if there is only one subform # sum the vectors Fj_petsc[j].zeroEntries()