Namespace issue with multiple adv/diff solvers

joss-paper/paper.md

@@ -72,20 +72,20 @@ The problems in global planetary dynamics and tectonics that `underworld3` is de
 
 ## Mathematical Framework
 
-`PETSc` provides a template form for the automatic generation of weak forms [see @knepley.etal.Achieving.2013]. The strong-form of the problem is defined through the functional $\cal{F}$ that expresses the balance between fluxes and unknowns:
+`PETSc` provides a template form for the automatic generation of weak forms [see @knepley.etal.Achieving.2013]. The strong-form of the problem is defined through the functional $\mathcal{F}$ that expresses the balance between fluxes, forces, and unknowns:
 
 \begin{equation}\label{eq:petsc-strong-form}
-\mathbf{F}(u) \sim \nabla \cdot F(u, \nabla u) - f(u, \nabla u) = 0
+\mathcal{F}(u) \sim \nabla \cdot F(u, \nabla u) - f(u, \nabla u) = 0
 \end{equation}
 
 The discrete weak form and its Jacobian derivative can be expressed as follows
 
 \begin{equation}\label{eq:petsc-weak-form}
- \mathbf{F}(u) \sim \sum_e \epsilon_e^T \left[ B^T W f(u^q, \nabla u^q) + \sum_k D_k^T W F^k (u^q, \nabla u^q) \right] = 0
+ \mathcal{F}(u) \sim \sum_e \epsilon_e^T \left[ B^T W f(u^q, \nabla u^q) + \sum_k D_k^T W F^k (u^q, \nabla u^q) \right] = 0
 \end{equation}
 
 \begin{equation}\label{eq:petsc-jacobian}
- \mathbf{F}'(u) \sim \sum _e \epsilon _{e^T}
+ \mathcal{F}'(u) \sim \sum _e \epsilon _{e^T}
                 \left[ \begin{array}{cc}
                     B^T  & D^T \\
                 \end{array} \right]
@@ -107,14 +107,14 @@ The discrete weak form and its Jacobian derivative can be expressed as follows
 The symbolic representation of the strong-form that is encoded in `underworld3` is:
 
 \begin{equation}\label{eq:sympy-strong-form}
-\underbrace{ \Bigl[ {D u}/{D t} \Bigr]}_{\dot{f} }
+\underbrace{ \Bigl[ {D u}/{D t} \Bigr]}_{\dot{u} }
 -\nabla \cdot \underbrace{\Bigl[ \mathrm{F}(u, \nabla u) \Bigr]}_{\mathbf{F}}
--\underbrace{\Bigl[ \mathrm{H}(\mathbf{x},t) \Bigr]}_{f}
+-\underbrace{\Bigl[ \mathrm{H}(\mathbf{x},t) \Bigr]}_{\mathbf{h}}
 = 0
 \end{equation}
 
 This symbolic form (\ref{eq:sympy-strong-form})
-contains material / time derivatives of the unknowns which are not present in the `PETSc` template because, after discretisation, these simplify to produce terms that are combinations of fluxes and flux history terms (which modify $F$) and forces (which modify $f$). In `underworld3`, the user interacts with the time derivatives themselves and `sympy` combines all the flux-like terms and all the force-like terms just prior to forming the Jacobians and compiling the `C` functions.
+contains material / time derivatives of the unknowns which are not present in the `PETSc` template because, after discretisation, these simplify to produce terms that are combinations of fluxes and flux history terms (which modify $F$) and forces (which modify $h$). In `underworld3`, the user interacts with the time derivatives themselves and `sympy` combines all the flux-like terms and all the force-like terms just prior to forming the Jacobians and compiling the `C` functions.
 
 # Discussion
 

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -25,10 +25,8 @@ class SolverBaseClass(uw_object):
     This class is not intended to be used directly
     """
 
-
     def __init__(self, mesh):
 
-
         super().__init__()
 
         self.mesh = mesh
@@ -669,8 +667,11 @@ class SNES_Scalar(SolverBaseClass):
         # f0 = sympy.Array(self._f0).reshape(1).as_immutable()
         # F1 = sympy.Array(self._f1).reshape(dim).as_immutable()
 
-        f0  = sympy.Array(uw.function.fn_substitute_expressions(self.F0.sym)).reshape(1).as_immutable()
-        F1  = sympy.Array(uw.function.fn_substitute_expressions(self.F1.sym)).reshape(dim).as_immutable()
+        # f0  = sympy.Array(uw.function.fn_substitute_expressions(self.F0.sym)).reshape(1).as_immutable()
+        # F1  = sympy.Array(uw.function.fn_substitute_expressions(self.F1.sym)).reshape(dim).as_immutable()
+
+        f0  = sympy.Array(uw.function.expressions.unwrap(self.F0.sym, keep_constants=False, return_self=False)).reshape(1).as_immutable()
+        F1  = sympy.Array(uw.function.expressions.unwrap(self.F1.sym, keep_constants=False, return_self=False)).reshape(dim).as_immutable()
 
         self._u_f0 = f0
         self._u_F1 = F1
@@ -1251,8 +1252,12 @@ class SNES_Vector(SolverBaseClass):
         # f0 = sympy.Array(self._f0).reshape(1).as_immutable()
         # F1 = sympy.Array(self._f1).reshape(dim).as_immutable()
 
-        f0  = sympy.Array(uw.function.fn_substitute_expressions(self.F0.sym)).reshape(dim).as_immutable()
-        F1  = sympy.Array(uw.function.fn_substitute_expressions(self.F1.sym)).reshape(dim,dim).as_immutable()
+        # f0  = sympy.Array(uw.function.fn_substitute_expressions(self.F0.sym)).reshape(dim).as_immutable()
+        # F1  = sympy.Array(uw.function.fn_substitute_expressions(self.F1.sym)).reshape(dim,dim).as_immutable()
+
+        f0  = sympy.Array(uw.function.expressions.unwrap(self.F0.sym, keep_constants=False, return_self=False)).reshape(dim).as_immutable()
+        F1  = sympy.Array(uw.function.expressions.unwrap(self.F1.sym, keep_constants=False, return_self=False)).reshape(dim,dim).as_immutable()
+
 
         self._u_f0 = f0
         self._u_F1 = F1
@@ -2060,10 +2065,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         ## and do these one by one as required by PETSc. However, at the moment, this
         ## is working .. so be careful !!
 
-
-        F0  = sympy.Array(uw.function.fn_substitute_expressions(self.F0.sym))
-        F1  = sympy.Array(uw.function.fn_substitute_expressions(self.F1.sym))
-        PF0 = sympy.Array(uw.function.fn_substitute_expressions(self.PF0.sym))
+        F0  = sympy.Array(uw.function.expressions.unwrap(self.F0.sym, keep_constants=False, return_self=False))
+        F1  = sympy.Array(uw.function.expressions.unwrap(self.F1.sym, keep_constants=False, return_self=False))
+        PF0  = sympy.Array(uw.function.expressions.unwrap(self.PF0.sym, keep_constants=False, return_self=False))
 
         # JIT compilation needs immutable, matrix input (not arrays)
         self._u_F0 = sympy.ImmutableDenseMatrix(F0)

src/underworld3/systems/ddt.py

@@ -38,7 +38,7 @@ def __init__(
         continuous: bool,
         swarm_degree: Optional[int] = None,
         swarm_continuous: Optional[bool] = None,
-        varsymbol: Optional[str] = r"u",
+        varsymbol: Optional[str] = None,
         verbose: Optional[bool] = False,
         bcs=[],
         order=1,
@@ -67,6 +67,9 @@ def __init__(
         else:
             self.swarm_continuous = swarm_continuous
 
+        if varsymbol is None:
+            varsymbol = rf"u_{{ [{self.instance_number}] }}"
+
         # meshVariables are required for:
         #
         # u(t) - evaluation of u_fn at the current time

src/underworld3/systems/solvers.py

@@ -1245,6 +1245,7 @@ def __init__(
                 vtype=uw.VarType.SCALAR,
                 degree=u_Field.degree,
                 continuous=u_Field.continuous,
+                varsymbol=u_Field.symbol,
                 verbose=verbose,
                 bcs=self.essential_bcs,
                 order=1,