underworldcode
UW constants
UW constants are sympy symbols that have an attached numerical value and a description. There is a subsitution operation attached to the class that will replace all constants with the numerical value and this is called in the function evaluation and JIT compiler. In future we might look to compile in PETSc constants in JIT.
The purpose of the constants is to make the "mathematically self describing" claim for uw3 closer to being true: we can write all parameters this way and keep all the help / view output symbolic instead of everything turning into crazy numbers.
I am not sure how this will interact with scaling.
"""
underworld `constants` are sympy symbols with attached
numeric values that are substituted into an underworld function
before evaluation. In sympy expressions, the symbol form is shown.
```{python}
alpha = uw_constant(
r'\\alpha',
value=3.0e-5,
description="thermal expansivity"
)
print(alpha.value)
print(alpha.description)
"""
If constants are nested in functions, sympy sees the substitution but if we also nest functions and constants, then we need to apply recursion as sympy can only see one layer deep in pulling out the atomic objects. This is not how I originally imagined the constants but this does mean we can define things like an effective viscosity with nested effective values (e.g. viscoelasticity with plasticity).
We could define a sub-expression class and have a sub class of constants which only allow float or int for their values
I've implemented the suggestion above - this functionality is now in the expressions class. Currently this does only work for single symbols (so we can't put in a symbolic tensor placeholder). Probably not that hard to do this, but the user-facing part might be a bit complicated to follow.
All tests are passing (perhaps we need more of them). The changes are to make the code better describe the mathematics that is is planning to solve and allow us to do this sort of thing:
The idea is also that we can try to make the equations more legible in the case where (for example) the viscosity is horribly non linear and full of sympy.piecewise branches
Introduce uw.function.derivative (should it live somewhere else to pair with Integration ?). Probably temporary - we should figure out the .diff from sympy but (see functions) this might be tricky to get 100% right.
Added various methods to identify if an expression has non-constant atoms anywhere in the chain. If not it can be considered a constant too and we just leave this alone as long as possible.
The sympy.Derivative machinery will happily make the Jacobians but they are nested matrices and we need to order them to match the petsc expectation. It may be easiest to flatten them, or apply the derivatives in the order we want.
Most recent commit addresses a number of quite complicated issues including some really tricky bugs that have been plaguing us when it comes to testing where multiple meshes are defined.
sympy does all sorts of magic behind the scenes on demand and it seems that it must do some caching that we cannot exactly control (i.e. sympy.core.cache.clear_cache() does not seem to be enough in the places I tried). The end result is that variables / expressions / functions that have the same name in sympy do not get replaced when they refer to a new variable - which occurs in the tests when we switch meshes. One workaround is to extend the name space in sympy to add some invisible characters to the symbol name for each mesh. I think this may be related to some of our memory / object management issues and needs more investigation.
A second caching issue arises from our mesh.evaluate code when we do not clear the cache (of data values) after accessing a mesh Variable. This has now been fixed.
@julesghub, @knepley - all the tests work and the current branch leans much more heavily on sympy than before. I realise we are doing more than @knepley's notion of the sympy functions as decorators or labels on our PETSc objects. It is almost the other way around because we are now really doing all the composing of operators in sympy itself.
You will see that I changed the python problem descriptions to simply be a collection of expressions to state F0 and F1 (corresponding to the relevant functions in PETSc) (e.g. for Stokes):
@property
def uf0(self):
DuDt = self.Unknowns.DuDt
f0 = uw.function.expression(
r"\mathbf{f}_0\left( \mathbf{u} \right)",
-self.bodyforce + self.rho * DuDt.bdf(1) / self.delta_t,
"NStokes pointwise force term: f_0(u)",
)
...
@property
def uF1(self):
dim = self.mesh.dim
DFDt = self.Unknowns.DFDt
if DFDt is not None:
# We can flag to only do this if the constitutive model has been updated
DFDt.psi_fn = self._constitutive_model.flux.T
F1 = uw.function.expression(
r"\mathbf{F}_1\left( \mathbf{u} \right)",
DFDt.adams_moulton_flux() - sympy.eye(self.mesh.dim) * (self.p.sym[0]),
"NStokes pointwise flux term: F_1(u)",
)
# Is the else condition useful - other than to prevent a crash ?
# Yes, because then it can just live on the Stokes solver ...
else:
...
the object uF0 is an uw.expression and can be differentiated with uw.function.derivative(). At the moment this is only partly implemented - Ideally it would return the Jacobians as uw.expressions so that we can see them both symbolically and as objects to pass to the compiler. This is a bit tedious but not complicated.
I also fixed up the viscoelastic / plastic constitutive model in the recent commits. This now uses the expressions to allow us to keep things like $\eta_\textrm{eff}$ as symbols in stokes.view() and other introspections.
uw.expressions know if they can be evaluated to pure constants (float, int or sympy.Number) and they are typically left alone in the printed symbols so that we don't have sympy going around combining them into long floating point numbers.
I have not tested the viscoelastic-plastic branch. There may be some requirements to express all of the max / min as approximations that are smooth for all the usual reasons.
@julesghub - we discussed today that we may want to harmonise the implementation of these expressions with the current implementation of uw.functions so that both use the sympy undefined function class as a starting point. This is more complicated but it might well be worth the effort. I suggest we merge the current code and work on how we can migrate the base class later.
Some updates on the NS benchmark case for @jcgraciosa -
-
Fixed a couple of bugs where the compilation was not getting properly updated for changes in timestep and density. This is a direct result of caching the expressions and I fixed this by expanding all the way down to the floating point values before hashing / caching.
-
I've been running this with P2/dP1 elements and the pressure field values seem much closer to the expected benchmark values than using the continuous pressures. Also, very sensitive to the resolution around the inclusion. Worth running again with this branch, I think.
Thanks, Louis. I'll re-test the models using this branch and the recommended settings.
Most recent commit: fixing the problem reported ages ago by Neng. When we deform the mesh, we had to hard-reset the solver. With this commit, the solver checks various cache-breaking conditions and resets itself accordingly. Tests passing.