diff --git a/joss-paper/paper.md b/joss-paper/paper.md index 00f58938..9c66ffd5 100644 --- a/joss-paper/paper.md +++ b/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 diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 8444c3b6..c4c9c3b6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/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) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 1e964293..ba1ec334 100644 --- a/src/underworld3/systems/ddt.py +++ b/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 diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index c48e0aab..c9ae12a2 100644 --- a/src/underworld3/systems/solvers.py +++ b/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,