Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Different ways of setting free-slip BCs have differing results #149

Closed
jcgraciosa opened this issue Jan 22, 2024 · 15 comments
Closed

Different ways of setting free-slip BCs have differing results #149

jcgraciosa opened this issue Jan 22, 2024 · 15 comments

Comments

@jcgraciosa
Copy link
Contributor

This is the issue I mentioned earlier in the meeting and posting here for documentation.

If I use different ways (shown below) of setting free-slip BCs then I get different results.

Method 1:

stokes.add_dirichlet_bc((0.0,), "Bottom", (1,))
stokes.add_dirichlet_bc((0.0,), "Top", (1,))
stokes.add_dirichlet_bc((0.0,), "Left", (0,))
stokes.add_dirichlet_bc((0.0,), "Right", (0,))

Method 2:

stokes.add_dirichlet_bc((sympy.oo,0.0), "Bottom")
stokes.add_dirichlet_bc((sympy.oo, 0.0), "Top")
stokes.add_dirichlet_bc((0.0,sympy.oo), "Left")
stokes.add_dirichlet_bc((0.0,sympy.oo), "Right")

Results:

Cartesian convection:
Resolution = 16:
Method 1:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -29.32204257096935
max. v (component): 23.151495282044586

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -27.873897927155742
max. v (component): 23.067035407567737

Resolution = 32:
Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINE_SEARCH iterations 8
min. v (component): -6.316221446999447
max. v (component): 4.949063245632038

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -27.820669550518556
max. v (component): 23.160314007131472

Resolution = 64
Method 1:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -29.398087782064863
max. v (component): 23.421398167185224

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -27.964331720326253
max. v (component): 22.900686663482414

SolC benchmark:
Resolution = 16:
Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
min. v (component): -2.738612787525831
max. v (component): 0.8660254037844388

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -0.023157758823695623
max. v (component): 0.02315529270672989

Resolution = 32:
Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
min. v (component): -2.738612787525831
max. v (component): 0.8660254037844388

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -0.023190411986986556
max. v (component): 0.02319011766071459

Resolution = 64:
Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
min. v (component): -2.738612787525831
max. v (component): 0.8660254037844388

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 2
min. v (component): -0.023198938179237818
max. v (component): 0.02319896279134057

Additional note:

  • used default petsc settings and same stokes solver tolerance value
@bknight1
Copy link
Member

The DIVERGED_LINEAR_SOLVE error in SolC suggests the BC isn't applied correctly.

if you try:

stokes.add_dirichlet_bc((0.0,0.0), "Bottom", (1,))
stokes.add_dirichlet_bc((0.0,0.0), "Top", (1,))
stokes.add_dirichlet_bc((0.0,0.0), "Left", (0,))
stokes.add_dirichlet_bc((0.0,0.0), "Right", (0,))

both method 1 and method 2 should produce the same result ( may have to be something like sympy.Matrix([0.,0.]) instead of (0.0,0.0) )

@lmoresi
Copy link
Member

lmoresi commented Jan 22, 2024 via email

@lmoresi
Copy link
Member

lmoresi commented Jan 23, 2024

We discussed switching from

stokes.add_dirichlet_bc((sympy.oo,0.0), "Bottom")

to this:

stokes.add_dirichlet_bc((None,0.0), "Bottom")

But the problem with this choice is that we can't pass in a sympy.Matrix with that form because

sympy.Matrix([None, 1, None]) is deprecated

1: SymPyDeprecationWarning:

non-Expr objects in a Matrix is deprecated. Matrix represents
a mathematical matrix. To represent a container of non-numeric
entities, Use a list of lists, TableForm, NumPy array, or some
other data structure instead.

See https://docs.sympy.org/latest/explanation/active-deprecations.html#deprecated-non-expr-in-matrix
for details.

This has been deprecated since SymPy version 1.9. It
will be removed in a future version of SymPy.

@julesghub @bknight1 @jcgraciosa - what do you think ?

@bknight1
Copy link
Member

Hmmm that's unfortunate but at least we've found out now.

I like sympy.oo under the hood for cancelling stuff out but not a big fan of it being user-facing. To me it's not intuitive that you'd use sympy.oo to have that component unconstrained. I think something like method 1 above is more intuitive, if the component is not set then it's unconstrained. I don't think it's that big of a deal though as once you've seen an example you'd know how to set free slip/no slip boundaries etc

@lmoresi
Copy link
Member

lmoresi commented Jan 23, 2024

I think it would be fine to use the UWGEO style bc description [None, 0, None] and sanitise it for sympy once we receive it in the solver. I like the idea that we are careful to specify the whole boundary in one go as it is very easy to specify inconsistent bc's otherwise.

@bknight1
Copy link
Member

So convert it to a sympy matrix under the hood, and covert None to sympy.oo, or something along those lines? I think that also sounds good!

@lmoresi
Copy link
Member

lmoresi commented Jan 24, 2024

Yes, basically. We currently do some processing to

  1. make sure it is a sympy.Matrix type
  2. determine which components are sympy.oo to construct a components list.

We do ultimately need a matrix for compilation as that's now assumed in the JIT code. We could leave it as sympy.oo because the compiler puts something useful into the C code when it sees that and it might help for debugging.

The real question is how to provide a pretty-printed view for the users. Perhaps we use a sympy symbol with a display form that is the empty set.

@NengLu
Copy link
Contributor

NengLu commented Jan 25, 2024

We may also need to include an extended boundary condition setting to make it possible to add extra node index sets (e.g., some internal nodes) flexibly.

In UWGeo style, it is like:
set_velocityBCs(left=[0,None], right=[0,None] top=[None,0], bottom=[None,0],nodeSets=[[fn,0],interface_nodes])

Now, it could be done in uw3 by embedding the mesh with an internal surface which is added as a label by the gmsh.
stokes.add_natural_bc((0,0), "Internal",(0,))

Could we make it more direct in petsc_generic_snes_solvers. _setup_discretisation() to add nodes in the boundary condition beside the "label"?

@lmoresi
Copy link
Member

lmoresi commented Jan 25, 2024

Have a look at the Kernel examples which do use an internal boundary.

Note that the PETSc interface for bc's is a bit hampered by its history and natural v essential bcs are a little bit different. At the moment we specify the natural bc's through a row matrix for all components because it is not entirely obvious what the components mean (to the user) on these embedded surfaces.

That's why I want to move to the same approach for other bc's.

@lmoresi
Copy link
Member

lmoresi commented Feb 21, 2024

@jcgraciosa - Issue noted, discussion above, can you finalise the proposed interface changes and we'll decide how to implement them.

@jcgraciosa
Copy link
Contributor Author

As a user, I think it's more intuitive to do the (None, 0., None) style. Also, since this is how BCs are set in UWGeo, sticking to it would be easier for long-time users (or those who need to use both uw3 and UWGeo) to transition.

@julesghub
Copy link
Member

Yep - please go with the UWGeo style @jcgraciosa

@lmoresi
Copy link
Member

lmoresi commented Feb 23, 2024

Internally, the sympy Matrix representation will need to use either sympy.oo or sympy.nan because None is not permitted as an entry in a sympy Matrix. I suggest a conversion to sympy.oo because it prints out better if you ask what boundary conditions have been set.

@jcgraciosa jcgraciosa removed their assignment Feb 29, 2024
jcgraciosa pushed a commit to jcgraciosa/underworld3 that referenced this issue Mar 6, 2024
@lmoresi
Copy link
Member

lmoresi commented Apr 3, 2024

@jcgraciosa - close if completed, please !

@lmoresi lmoresi closed this as completed Apr 3, 2024
@lmoresi lmoresi reopened this Apr 3, 2024
@lmoresi
Copy link
Member

lmoresi commented Apr 3, 2024

Sorry - not mine to close.

julesghub pushed a commit that referenced this issue Apr 24, 2024
julesghub added a commit that referenced this issue Jun 4, 2024
* Weird scale-dependenct swarm-location bug

* Add ._utils file

File to place functions that may be useful. First function is to gather data on the root proc which could also be broadcast to all processors.

* Adding a solver base class to simplify the documentation. Adding class method for ipython help so that we can ask about solvers before using them.

* Tests passing again

* Lots of changes needed to get the viscoelastic models close to working. These are in place now and I am going back to try to get the tests and notebooks working again.

* Ooops, some python 3.10 stuff that is probably overly cutting edge

* Create docker-image.yml

Add a action to build and push a docker to dockerhub

* Update docker-image.yml

invalid tag name fix

* Update docker-image.yml

* add branch name

* Update docker-image.yml

* Update docker-image.yml

Setting docker build on VEP branch only

* fix adv_dif solver

* Update build_uw_test.yml

Fix cython version to <3.0

* Fix divide by 0 error

Re-write estimate_dt to remove divide by 0 error

* Fix divide by 0 error

Re-write estimate_dt to remove divide by 0 error

* switch to pdoc (not pdoc3) which seems better supported and the math mode matches jupyter

* Preliminary SLCN implementation for Naver Stokes and generalising this to be a D/DT module.

* Bugz introduced in solvers.py

* tiny bug in cython code

* Needed to change callback casts for new Cython

* restoring the petsc swarm populate routines

* Updates for Cython 3 and latest PETSc

* Whitespace

* UnderworldFunction should not hold the mesh

* MeshVariable symbols should not hold the mesh

* MeshVariable should not hold the PetscFE

* Fix indentation mistake

* Respect argument to mesh.access()

* Do not create DS until all mesh variables are defined

* Added assert()

* Some geom mg changes

* Revert discretization changes until rewrite

* Rvert change for leaks until redesign

* Revert "MeshVariable symbols should not hold the mesh"

This reverts commit f70f898.

* I thought I might need DMSetLocalSection()

* Fix SaddlePoint field names

* Do not recreate the DM in the SaddlePoint consrtuctor

* Do not recreate the DS in the SaddlePoint constructor

* Must copy the discretization to coarse grids

* WARNING Turning off auxiliary vector

* SolCx works with additional fields defined / set (not used)

* Fix up MKMG - backwards compatible etc etc. Tests pass, bugs may just be hiding a bit better than before

* Quite a few post-knepley updates - see pull request for details

* Issues with the Stokes solver but Poisson is working correctly again.

* Placeholder commit - working nbc, ebc ... tests need some work

* Changes to get solvers working

Few changes to keep up with Louis changes. Working: AdvDiff_swarm, Darcy, Stokes.

Not Working:
- *Adv_diff*:

-> 1205 self._F_star_projection_solver.solve()
   1207 with self._nswarm_F.access(nX0):
   1208     nX0.data[...] = self._nswarm_F.data[...]
ValueError: could not broadcast input array from shape (33282,) into shape (73728,)

Notes:
- Had to update BC in workflow, issue with int64 (python) and long (petsc) np array on mac. Had to set np array to int32.
- Have to set each BC individually

* Fix up semi-Lagrange manager and take a small step towards eradicating the bat

* Some fixes for anti-batman strategy, adv-diff test, swarms clone cell dm. More testing required esp in parallel

* Updating batman / semiLagrange fixes

* Couple of updates to get NavierStokesSLCN running

* Remove pycache directories too in clean.sh

* Some memory cleaning.
Underworldfunction no longer hold a mesh. Instead they reference the
mesh via the MeshVariable associated with the uw_fn.

* Checking if petsc4py is installed under the PETSC_DIR directory.
If not a warning is produced - good enough for now.

* evalf() function need to use the meshvariable.mesh too, as per previous
commit

* For add_dirichlet_bc() fix 'component' for 'components'.

* Remove some obsolete lines from solves.py

* Checking in - passing tests but incomplete

* Fixing some maths-in-comments bugs (r-strings)

* WIP ideas only, delete patch files please

* Update solvers.py

Add in CM to Stokes (NS inherits from Stokes) and SNES_Poisson (Adv_diff inherits from Poisson)

Should be backwards compatible in current state (🤞)

TODO:
* Update constitutive_models base class to include: DfDt and DuDt vars to be passed from the solver setup to the CM
* multimaterial bits

* Merge commit 1

* Updates to CM

Pass through DuDt & DFDt to CM for use.
Now attaches unknown to the CM under the hood

* Batman returns !!!!

* Solver & CM updates

Trying to make sure DFDt and DuDt are the same in the solver and CM for updating values that are used in both

* Update solvers.py

Messed up the order of variables when passing them to the CM.

* CM & solver stuff AGAIN

Changes to use the lazy eval approach.

* Update petsc_compat layer to remove warnings - would be nice to automat this process from petscds.h

* Update solvers to remove class inheritance

passing tests....
working from UW3-benchmarks:
- Darcy (without g term)
- Stokes
- NS - SLCN
- AdvDiff-SLCN
-AdvDiff-swarm
Not working:
- NS-swarm (solver crashes)

* Prototype of the Unknowns sub-class for discussion purposes.

* closer but one test failing

* Semi Lagrange / fully Lagrangian history terms now have a consistent interface. There is an issue with the sequencing that I haven't figured out yet so the full Lagrangian version is not updating correctly.

* removing obsolete duplicate files

* Minor interface tweet for Adv Diff solver

* Renaming the petsc options so linux+petsc find's the options

* fixing a batman bug

* Moving unknowns into the Unknowns class: u, L=grad(u), E=1/2(L + L^T), W=1/2(L-L^T)

* Prep for merge

* tidy up

* Version number

* Fixing up constitutive models and solvers so tests pass when using the solver Unknowns class ... VE Stokes  needs work.

* VE Stokes as a separate solver type - I give up on trying to make a clever way to add VE history tracking. This is OK. tests added and passing. I won't vouch for Navier Stokes yet. That needs refactoring in line with the constitutive model updates.

* quick syntax fix for Examples-Convection

* validate_solver() is back, removing _Unknowns from consitutive model as we have getters and setters for this

* Couple of changes needed to get Navier-Stokes (SLCN) over the line with all the other changes we did. This is not really well tested.

* Add some tests

AdvDiff & Darcy tests in cartesian

* Change for petsc main compatibility projectCoordinates() -> setCoordinateDisc()

with argument setting setCoordinateDisc(..., project=False) - is quicker i understand.

removing obsolete petsc_dm_project_coordinates() from
petsc_discretisation.pyx

* Read MSH file too

* Fixing solver and visualization calls

* Fixing solver and visualization calls

* Fixing solver and visualization calls

* Update for underworld3-0.9b conda package

* conda improvements for linux

* Fix for platform agnostic gxx compiler

* set conda to build from development branch for now

* Debugging the natural BC / surface integral terms

* Fix up tests (down-size the AdvDiff for now). Some merges required.

* from mpi4py import MPI before petsc

Solves issue #140

Cleaning up some unused python imports too.

* Surface integrals working with clobbered version of petsc for free-slip bcs using penalty approach. Example uploaded in WIP as a demo of this and to show what happens with current petsc main branch. Running version of cylindrical/annulus convection model with free-slip upper bc.

* Changes to SetSNESLocalFEM call pattern ...

* Ooops - stats() is collective !

* Getting the internal-heated discworld to run with Free-Slip (flux limiting on the advection scheme with a non-linear diffusivity)

* Move the surface integral example out of WIP

* Using PETSc normals via a new vector mesh.Gamma which maps to petsc_n in surface integrals. See Ex_Stokes_Ellipse for an example

* Tidying up - adding the channel flow / irregular base model

* Setting default inner radius for spherical meshes to be Earth like (0.547 not 0.3). Adding a regional spherical mesh since we can now do the side bcs easily. Should consider a bucket mesh too, and one with cuts along any great circle.

* Fixing Darcy Solver and updating example. Removing some old notebooks.

* Adding an example that uses the RegionalCap Mesh. Why not since we have the impermeable boundary thing working well.

* Fix broken pointwise function test - see comments

* Update meshing.py

Fix for #147

* Trying to recreate kernel models in 2D / annulus

* Fix up broken tests and fix Darcy solver

* Adding Cartesian Kernel prototype example (strange mesh, sorry !) - starting point for Parsons & Daly benchmark

* Updating the Parsons & Daly Kernels model. Not benchmarked, just a demo.

* Updates to allow petsc options from cmd line. Disc Convection model and visualisation. Change snes_conergence test checking (now responsibility of the solver)

* Checking a few things for parallel surface integrals

* Simplify boundary condition logic by relying on mesh.boundary Enum for values

* Forgot to update Stokes

* Adding a couple of workaround meshes to carry on benchmarking while we sort out the boundary integral stuff.

* PETSc.FE().createDefault(..., PETSc.COMM_SELF) - thanks Matt

* Found a bug in advection swarm

* Moving the surface integrals example as this is now in PETSc main. Clean up some things.

* Adjusting SolCx to compare boundary condition types

* Null space calculation annulus

* Trying to fix the solver-then-new-var batman-returns bug. Update to some new petsc option names. Examples with projected normals for irregular boundaries.

* Fixing Darcy Notebook. Script to run tests independently. Adv-Diff interpolation still problematic

* Changes to src structure (re-built from jg branch)

* Union - cython error in build system

* cython is more fussy in this environment

* Out of date cython & petsc

* Petsc minor release API change ... honestly WTF

* removing obsolete line

* Putting import MPI first (again) and improving swarm pytest with a
fixture

* lm src style (#159)

* Changes to src structure (re-built from jg branch)

* Union - cython error in build system

* cython is more fussy in this environment

* Out of date cython & petsc

* Petsc minor release API change ... honestly WTF

* removing obsolete line

* Putting import MPI first (again) and improving swarm pytest with a
fixture

---------

Co-authored-by: Julian Giordani <house.of.jules@gmail.com>

* compatibility thing ... trying to get the build to work on old petsc. But this will probably not pass tests because it's the old petsc

* Fix calling pattern

* Update discretisation.py

* Drop tests that we know will fail with the conda-forge version of petsc

* Re-ordering tests by category

* Compatible petsc4py is available on conda-forge

* Update docs to find the new path (duh !)

* Remember to actually install a recent pdoc

* Update deployment path

* Some small re-arrangements to make the pdoc API docs look a bit cleaner.

* Take out the "numpy" from sympy lambdify - this automatically includes scipy if available and allows compilation of special functions.

* Inconsistencies in evalf / evaluate in function.pyx. Fixing that and also allowing scalar / vector functions again. Tests passing.

* Visualisation module needs small update to match

* Checking in the 1d block analytic solution for the Adv Diff problem in 2D (really 1D).

* Updated Example notebook for AdvDiff

* Update issue templates (#162)

First cut at a bug report. Done #102 ✅

* Updating the 1D advection example with a diagram

* Found some bugs in the way we determine element sizes. Made this a kd-tree thing to step away from the PETSc FVM code. Fix up some remaining NS issues. Moved the DDt objects into systems where they belong (not in swarm.py). Tidied that lot up and fixed a couple of tests.

* update petsc 3.20.5 for build / test

* Not petsc4py, apparently

* Update micromamba

* This micromamba is on a different path

* Still have this compat / wrapper issue

* mpmath version error (good grief)

* openGL errors in gmsh

* workflow / gmsh

* Updates for current petsc-main

* Quick Fix for proxy variable interpolation on sparse swarms

* Yeah, OK, nnn is the limit. Not  permanent fix.

* Test binder (#171)

* Yeah, OK, nnn is the limit. Not  permanent fix.

* Add binder files

* uw3_env fails (use std petsc for now)

* Full install into binder python

* Use the petsc4py paths to PETSc

* gmsh is fussy on binder (obviously)

* going the other way - add everything gmsh MIGHT need

* Fixing meshVariable.clone - not properly initialised. Adding simplify to visualisation evaluation.

* Adding the simplify to all function evaluations (helps with matrix * piecewise for example, or matrix components)

* try adding a start file

* try conda activate script from postbuild

* typo / conda

* uw3_env - 0.0

* typo - I forgot to add uw3_env

* move all the uw3 stuff to postbuild

* postBuild - update

* uw3_env 0.1

* uw3_env 0.2

* (readme)

* uw3_env 0.3

* uw3_env 0.4

* uw3_env 0.5

* uw3_env 0.6

* uw3_env 0.7

* ue3_env 0.8

* uw3_env 0.8

* uw3_env 0.9

* backing off - uw3_env

* fails for py3.10 but not 3.11 ?

* py 3.11 builds correctly.. Now debug environment

* more info

* conda is not activated in postBuild

* and ...

* simplify scripts - binder is now working

* accidentally deleted the uw build - put it back

* Some changes for the merge dev -> main

* moving some docs around, switch to quarto, fix mesh self-documentation

* Bug in mesh.view

* Bc stacker (#177)

* Stacking bcs in a single label to simplify internal implementation. Tests pass etc

* Careful with empty labels

* Gadjets for manipulating labels as masks

* Adding All_Boundaries to bc enum so UW can display it

* small fix to bc masks

* verbose - bc problems

* empty label problem

* Check that differently

* verbose snes solve

* verbose NS solver

* evalf option for ddt

* Add some verbose options to ddt

* also in advection

* Pulling some of those calls back out

* Remove some print calls

* Allow ns.solve to over-ride self.order which allows

* typo

* update NS solver

* revert

* nswarm changes

* re-introduce timestep averaging for ddt / semi-lagrange

* propagate the physical dt

* adding the dt_physical as a NS property

* bc_mask for NS Solver / advection term

* push through changes to NS

* sympy.matrix bugs

* another sympy typon

* sympy typo 3

* sympy 4

* sympy5

* bc_mask - put this in the nodal swarm so boundary points can be left fixed even for non-zero vbcs

* ddt - typo

* access required in swarm

* np shape issue

* parallel issues - nswarm

* oops - shape again

* sympy function issue

* verify function

* evalf test

* doing the masks in numpy not sympy

* Particle advection with sub-stepping (nswarm only)

* h5 viewer for mesh

* clean up, clean up, everybody everywhere

* Test JIT directly (not via stokes setup)

* substepping for regular swarm

* cache / JIT arguments

* Make test_0004 clean up better (before and after)

* clean up swarm code

* binder: use 3.21.0 from conda-forge

* .binder update

* Try not pinning python

* Fix for issue #179 (#180)

issue #179

Co-authored-by: Juan Carlos Graciosa <juan.graciosa@monash.edu>

* move all documentation / jupyterbook to its own repo

* Park the doc folder as well

* Passive swarm fixes apr24 (#183)

* Fixes for swarm substepping

* Checking

* Monitor 1

* breaking in numpy - check that instead

* typo

* Odd

* interaction between swarm.access and swarm.dm.getLocalSize()

* IndexError also possible ...

* typo - fix

* fix inner/outer surface gmsh issue (#184)

* Ugh - I added a huge, hairy bug to the advection solver (ddt.py) which probably broke all the NS and AD benchmarks. Suggest we add a test for this based on the 1D adv-diffusion solver example which finally picked up the problem.

* Fix tests - adv diffusion

* typo in swarm.py that breaks parallel advection

* Marking merged versions of SLCN / Lagrangian ADV-DIFF solvers as deleted (they are currently not imported at all).

* Fix pyvista issue while saving a screenshot (#188)

changed jupyter_backend to trame

* Added function to handle None inputs in bcs

* Issue #149 fix

* Only handle iterable BC inputs

* Simplifying BC conversion code.

All conditions are now turned into a list before conversion.
Modifying test to test the conversion code.

* conda recipe fixes (#190)

conda recipe updates

* Updating the interpolation routines. (#191)

Updating the interpolation routines.

These c functions are direct copies of the petsc routines
DMInterpolation_SetUp and DMInterpolation_Evaluate, only
with owning_cell information provided and used. Until petsc can
optionally use 'owning_cells' we will have to manually update the c code
like this commit.

* Julesghub/mpi update (#192)

* importing underworld3.mpi earlier on to avoid circular imports

* Only get the communicator via underworld3.mpi.comm

* fix pyvista issue in binder (#194)

fix pyvista issue in binder, check whether it is running in binder or not

---------

Co-authored-by: Louis Moresi <louis.moresi@anu.edu.au>
Co-authored-by: Ben Knight <55677727+bknight1@users.noreply.github.com>
Co-authored-by: Matthew G. Knepley <knepley@gmail.com>
Co-authored-by: Tyagi <gthyagi@gmail.com>
Co-authored-by: jcgraciosa <12045001+jcgraciosa@users.noreply.github.com>
Co-authored-by: Juan Carlos Graciosa <juan.graciosa@monash.edu>
Co-authored-by: Thyagarajulu Gollapalli <44185316+gthyagi@users.noreply.github.com>
julesghub added a commit that referenced this issue Jun 4, 2024
* Adding a solver base class to simplify the documentation. Adding class method for ipython help so that we can ask about solvers before using them.

* Tests passing again

* Lots of changes needed to get the viscoelastic models close to working. These are in place now and I am going back to try to get the tests and notebooks working again.

* Ooops, some python 3.10 stuff that is probably overly cutting edge

* Create docker-image.yml

Add a action to build and push a docker to dockerhub

* Update docker-image.yml

invalid tag name fix

* Update docker-image.yml

* add branch name

* Update docker-image.yml

* Update docker-image.yml

Setting docker build on VEP branch only

* fix adv_dif solver

* Update build_uw_test.yml

Fix cython version to <3.0

* Fix divide by 0 error

Re-write estimate_dt to remove divide by 0 error

* Fix divide by 0 error

Re-write estimate_dt to remove divide by 0 error

* switch to pdoc (not pdoc3) which seems better supported and the math mode matches jupyter

* Preliminary SLCN implementation for Naver Stokes and generalising this to be a D/DT module.

* Bugz introduced in solvers.py

* tiny bug in cython code

* Needed to change callback casts for new Cython

* restoring the petsc swarm populate routines

* Updates for Cython 3 and latest PETSc

* Whitespace

* UnderworldFunction should not hold the mesh

* MeshVariable symbols should not hold the mesh

* MeshVariable should not hold the PetscFE

* Fix indentation mistake

* Respect argument to mesh.access()

* Do not create DS until all mesh variables are defined

* Added assert()

* Some geom mg changes

* Revert discretization changes until rewrite

* Rvert change for leaks until redesign

* Revert "MeshVariable symbols should not hold the mesh"

This reverts commit f70f898.

* I thought I might need DMSetLocalSection()

* Fix SaddlePoint field names

* Do not recreate the DM in the SaddlePoint consrtuctor

* Do not recreate the DS in the SaddlePoint constructor

* Must copy the discretization to coarse grids

* WARNING Turning off auxiliary vector

* SolCx works with additional fields defined / set (not used)

* Fix up MKMG - backwards compatible etc etc. Tests pass, bugs may just be hiding a bit better than before

* Quite a few post-knepley updates - see pull request for details

* Issues with the Stokes solver but Poisson is working correctly again.

* Placeholder commit - working nbc, ebc ... tests need some work

* Changes to get solvers working

Few changes to keep up with Louis changes. Working: AdvDiff_swarm, Darcy, Stokes.

Not Working:
- *Adv_diff*:

-> 1205 self._F_star_projection_solver.solve()
   1207 with self._nswarm_F.access(nX0):
   1208     nX0.data[...] = self._nswarm_F.data[...]
ValueError: could not broadcast input array from shape (33282,) into shape (73728,)

Notes:
- Had to update BC in workflow, issue with int64 (python) and long (petsc) np array on mac. Had to set np array to int32.
- Have to set each BC individually

* Fix up semi-Lagrange manager and take a small step towards eradicating the bat

* Some fixes for anti-batman strategy, adv-diff test, swarms clone cell dm. More testing required esp in parallel

* Updating batman / semiLagrange fixes

* Couple of updates to get NavierStokesSLCN running

* Remove pycache directories too in clean.sh

* Some memory cleaning.
Underworldfunction no longer hold a mesh. Instead they reference the
mesh via the MeshVariable associated with the uw_fn.

* Checking if petsc4py is installed under the PETSC_DIR directory.
If not a warning is produced - good enough for now.

* evalf() function need to use the meshvariable.mesh too, as per previous
commit

* For add_dirichlet_bc() fix 'component' for 'components'.

* Remove some obsolete lines from solves.py

* Checking in - passing tests but incomplete

* Fixing some maths-in-comments bugs (r-strings)

* WIP ideas only, delete patch files please

* Update solvers.py

Add in CM to Stokes (NS inherits from Stokes) and SNES_Poisson (Adv_diff inherits from Poisson)

Should be backwards compatible in current state (🤞)

TODO:
* Update constitutive_models base class to include: DfDt and DuDt vars to be passed from the solver setup to the CM
* multimaterial bits

* Merge commit 1

* Updates to CM

Pass through DuDt & DFDt to CM for use.
Now attaches unknown to the CM under the hood

* Batman returns !!!!

* Solver & CM updates

Trying to make sure DFDt and DuDt are the same in the solver and CM for updating values that are used in both

* Update solvers.py

Messed up the order of variables when passing them to the CM.

* CM & solver stuff AGAIN

Changes to use the lazy eval approach.

* Update petsc_compat layer to remove warnings - would be nice to automat this process from petscds.h

* Update solvers to remove class inheritance

passing tests....
working from UW3-benchmarks:
- Darcy (without g term)
- Stokes
- NS - SLCN
- AdvDiff-SLCN
-AdvDiff-swarm
Not working:
- NS-swarm (solver crashes)

* Prototype of the Unknowns sub-class for discussion purposes.

* closer but one test failing

* Semi Lagrange / fully Lagrangian history terms now have a consistent interface. There is an issue with the sequencing that I haven't figured out yet so the full Lagrangian version is not updating correctly.

* removing obsolete duplicate files

* Minor interface tweet for Adv Diff solver

* Renaming the petsc options so linux+petsc find's the options

* fixing a batman bug

* Moving unknowns into the Unknowns class: u, L=grad(u), E=1/2(L + L^T), W=1/2(L-L^T)

* Prep for merge

* tidy up

* Version number

* Fixing up constitutive models and solvers so tests pass when using the solver Unknowns class ... VE Stokes  needs work.

* VE Stokes as a separate solver type - I give up on trying to make a clever way to add VE history tracking. This is OK. tests added and passing. I won't vouch for Navier Stokes yet. That needs refactoring in line with the constitutive model updates.

* quick syntax fix for Examples-Convection

* validate_solver() is back, removing _Unknowns from consitutive model as we have getters and setters for this

* Couple of changes needed to get Navier-Stokes (SLCN) over the line with all the other changes we did. This is not really well tested.

* Add some tests

AdvDiff & Darcy tests in cartesian

* Change for petsc main compatibility projectCoordinates() -> setCoordinateDisc()

with argument setting setCoordinateDisc(..., project=False) - is quicker i understand.

removing obsolete petsc_dm_project_coordinates() from
petsc_discretisation.pyx

* Read MSH file too

* Fixing solver and visualization calls

* Fixing solver and visualization calls

* Fixing solver and visualization calls

* Update for underworld3-0.9b conda package

* conda improvements for linux

* Fix for platform agnostic gxx compiler

* set conda to build from development branch for now

* Debugging the natural BC / surface integral terms

* Fix up tests (down-size the AdvDiff for now). Some merges required.

* from mpi4py import MPI before petsc

Solves issue #140

Cleaning up some unused python imports too.

* Surface integrals working with clobbered version of petsc for free-slip bcs using penalty approach. Example uploaded in WIP as a demo of this and to show what happens with current petsc main branch. Running version of cylindrical/annulus convection model with free-slip upper bc.

* Changes to SetSNESLocalFEM call pattern ...

* Ooops - stats() is collective !

* Getting the internal-heated discworld to run with Free-Slip (flux limiting on the advection scheme with a non-linear diffusivity)

* Move the surface integral example out of WIP

* Using PETSc normals via a new vector mesh.Gamma which maps to petsc_n in surface integrals. See Ex_Stokes_Ellipse for an example

* Tidying up - adding the channel flow / irregular base model

* Setting default inner radius for spherical meshes to be Earth like (0.547 not 0.3). Adding a regional spherical mesh since we can now do the side bcs easily. Should consider a bucket mesh too, and one with cuts along any great circle.

* Fixing Darcy Solver and updating example. Removing some old notebooks.

* Adding an example that uses the RegionalCap Mesh. Why not since we have the impermeable boundary thing working well.

* Fix broken pointwise function test - see comments

* Update meshing.py

Fix for #147

* Trying to recreate kernel models in 2D / annulus

* Fix up broken tests and fix Darcy solver

* Adding Cartesian Kernel prototype example (strange mesh, sorry !) - starting point for Parsons & Daly benchmark

* Updating the Parsons & Daly Kernels model. Not benchmarked, just a demo.

* Updates to allow petsc options from cmd line. Disc Convection model and visualisation. Change snes_conergence test checking (now responsibility of the solver)

* Checking a few things for parallel surface integrals

* Simplify boundary condition logic by relying on mesh.boundary Enum for values

* Forgot to update Stokes

* Adding a couple of workaround meshes to carry on benchmarking while we sort out the boundary integral stuff.

* PETSc.FE().createDefault(..., PETSc.COMM_SELF) - thanks Matt

* Found a bug in advection swarm

* Moving the surface integrals example as this is now in PETSc main. Clean up some things.

* Adjusting SolCx to compare boundary condition types

* Null space calculation annulus

* Trying to fix the solver-then-new-var batman-returns bug. Update to some new petsc option names. Examples with projected normals for irregular boundaries.

* Fixing Darcy Notebook. Script to run tests independently. Adv-Diff interpolation still problematic

* Changes to src structure (re-built from jg branch)

* Union - cython error in build system

* cython is more fussy in this environment

* Out of date cython & petsc

* Petsc minor release API change ... honestly WTF

* removing obsolete line

* Putting import MPI first (again) and improving swarm pytest with a
fixture

* lm src style (#159)

* Changes to src structure (re-built from jg branch)

* Union - cython error in build system

* cython is more fussy in this environment

* Out of date cython & petsc

* Petsc minor release API change ... honestly WTF

* removing obsolete line

* Putting import MPI first (again) and improving swarm pytest with a
fixture

---------

Co-authored-by: Julian Giordani <house.of.jules@gmail.com>

* compatibility thing ... trying to get the build to work on old petsc. But this will probably not pass tests because it's the old petsc

* Fix calling pattern

* Update discretisation.py

* Drop tests that we know will fail with the conda-forge version of petsc

* Re-ordering tests by category

* Compatible petsc4py is available on conda-forge

* Update docs to find the new path (duh !)

* Remember to actually install a recent pdoc

* Update deployment path

* Some small re-arrangements to make the pdoc API docs look a bit cleaner.

* Take out the "numpy" from sympy lambdify - this automatically includes scipy if available and allows compilation of special functions.

* Inconsistencies in evalf / evaluate in function.pyx. Fixing that and also allowing scalar / vector functions again. Tests passing.

* Visualisation module needs small update to match

* Checking in the 1d block analytic solution for the Adv Diff problem in 2D (really 1D).

* Updated Example notebook for AdvDiff

* Update issue templates (#162)

First cut at a bug report. Done #102 ✅

* Updating the 1D advection example with a diagram

* Found some bugs in the way we determine element sizes. Made this a kd-tree thing to step away from the PETSc FVM code. Fix up some remaining NS issues. Moved the DDt objects into systems where they belong (not in swarm.py). Tidied that lot up and fixed a couple of tests.

* update petsc 3.20.5 for build / test

* Not petsc4py, apparently

* Update micromamba

* This micromamba is on a different path

* Still have this compat / wrapper issue

* mpmath version error (good grief)

* openGL errors in gmsh

* workflow / gmsh

* Updates for current petsc-main

* Quick Fix for proxy variable interpolation on sparse swarms

* Yeah, OK, nnn is the limit. Not  permanent fix.

* Test binder (#171)

* Yeah, OK, nnn is the limit. Not  permanent fix.

* Add binder files

* uw3_env fails (use std petsc for now)

* Full install into binder python

* Use the petsc4py paths to PETSc

* gmsh is fussy on binder (obviously)

* going the other way - add everything gmsh MIGHT need

* Fixing meshVariable.clone - not properly initialised. Adding simplify to visualisation evaluation.

* Adding the simplify to all function evaluations (helps with matrix * piecewise for example, or matrix components)

* try adding a start file

* try conda activate script from postbuild

* typo / conda

* uw3_env - 0.0

* typo - I forgot to add uw3_env

* move all the uw3 stuff to postbuild

* postBuild - update

* uw3_env 0.1

* uw3_env 0.2

* (readme)

* uw3_env 0.3

* uw3_env 0.4

* uw3_env 0.5

* uw3_env 0.6

* uw3_env 0.7

* ue3_env 0.8

* uw3_env 0.8

* uw3_env 0.9

* backing off - uw3_env

* fails for py3.10 but not 3.11 ?

* py 3.11 builds correctly.. Now debug environment

* more info

* conda is not activated in postBuild

* and ...

* simplify scripts - binder is now working

* accidentally deleted the uw build - put it back

* Some changes for the merge dev -> main

* moving some docs around, switch to quarto, fix mesh self-documentation

* Bug in mesh.view

* Bc stacker (#177)

* Stacking bcs in a single label to simplify internal implementation. Tests pass etc

* Careful with empty labels

* Gadjets for manipulating labels as masks

* Adding All_Boundaries to bc enum so UW can display it

* small fix to bc masks

* verbose - bc problems

* empty label problem

* Check that differently

* verbose snes solve

* verbose NS solver

* evalf option for ddt

* Add some verbose options to ddt

* also in advection

* Pulling some of those calls back out

* Remove some print calls

* Allow ns.solve to over-ride self.order which allows

* typo

* update NS solver

* revert

* nswarm changes

* re-introduce timestep averaging for ddt / semi-lagrange

* propagate the physical dt

* adding the dt_physical as a NS property

* bc_mask for NS Solver / advection term

* push through changes to NS

* sympy.matrix bugs

* another sympy typon

* sympy typo 3

* sympy 4

* sympy5

* bc_mask - put this in the nodal swarm so boundary points can be left fixed even for non-zero vbcs

* ddt - typo

* access required in swarm

* np shape issue

* parallel issues - nswarm

* oops - shape again

* sympy function issue

* verify function

* evalf test

* doing the masks in numpy not sympy

* Particle advection with sub-stepping (nswarm only)

* h5 viewer for mesh

* clean up, clean up, everybody everywhere

* Test JIT directly (not via stokes setup)

* substepping for regular swarm

* cache / JIT arguments

* Make test_0004 clean up better (before and after)

* clean up swarm code

* binder: use 3.21.0 from conda-forge

* .binder update

* Try not pinning python

* Fix for issue #179 (#180)

issue #179

Co-authored-by: Juan Carlos Graciosa <juan.graciosa@monash.edu>

* move all documentation / jupyterbook to its own repo

* Park the doc folder as well

* Passive swarm fixes apr24 (#183)

* Fixes for swarm substepping

* Checking

* Monitor 1

* breaking in numpy - check that instead

* typo

* Odd

* interaction between swarm.access and swarm.dm.getLocalSize()

* IndexError also possible ...

* typo - fix

* fix inner/outer surface gmsh issue (#184)

* Ugh - I added a huge, hairy bug to the advection solver (ddt.py) which probably broke all the NS and AD benchmarks. Suggest we add a test for this based on the 1D adv-diffusion solver example which finally picked up the problem.

* Fix tests - adv diffusion

* typo in swarm.py that breaks parallel advection

* Marking merged versions of SLCN / Lagrangian ADV-DIFF solvers as deleted (they are currently not imported at all).

* Fix pyvista issue while saving a screenshot (#188)

changed jupyter_backend to trame

* Added function to handle None inputs in bcs

* Issue #149 fix

* Only handle iterable BC inputs

* Simplifying BC conversion code.

All conditions are now turned into a list before conversion.
Modifying test to test the conversion code.

* conda recipe fixes (#190)

conda recipe updates

* Updating the interpolation routines. (#191)

Updating the interpolation routines.

These c functions are direct copies of the petsc routines
DMInterpolation_SetUp and DMInterpolation_Evaluate, only
with owning_cell information provided and used. Until petsc can
optionally use 'owning_cells' we will have to manually update the c code
like this commit.

* Julesghub/mpi update (#192)

* importing underworld3.mpi earlier on to avoid circular imports

* Only get the communicator via underworld3.mpi.comm

* fix pyvista issue in binder (#194)

fix pyvista issue in binder, check whether it is running in binder or not

* UW constants (#187)

* Add uw_constants functionality (available, but not covered in tests, plugged into jit compiler and function evaluation)

* uw_constants: make recursive so we can nest constants in other constants correctly. Also works for nests of functions, constants etc, so perhaps this is really a sub-expression class not a constant class. Fix typo (substitute)

* renaming to 'expressions' and providing capacity to substitute all or just one expression. ViscousFlow and ViscoPlastic are converted over. VEP not yet

* Converting the solvers so f0, F1 are supplied as properties and any side effects are handled there (e.g. loading RHS for associated projection solvers). This simplifies the flow and make f0, F1 available before solve time so these can be checked by the user beforehand. All tests passing (not all at once, unfortunately)

* Tidy up before merging

* Surface integral tests should go back in

* Fixing up Viscoelasticity using expressions. Need to add meaningful tests. No word on whether VEP produces proper Jacobians or if we need to use a smooth version of Min()/Max()

* Managing differentiation of quantities with expressions. Ideally we intercept sympy.diff and make substitutions for non-constant expressions first (sympy sees all the symbols as constants). Not trivial because of sympy's approach to cloning objects.

* Fixes for VEP - DDt / functions. Better type-checking in expressions. Update object viewer for meshvariables (so it is actually helpful)

* Updates to fix sympy and caching issues

* Fixes for problems in Navier-Stokes

* Caching mechanism for JIT needs to expand expressions to see if their content has been updated (esp for things like timestep)

* Not sure u.degree - 1 or just 1 for the interpolation degree of the flux field terms.

* Implement penalty term in NS (go inheritence, go !!)

* Fix read_timestep function for discontinuous variables.

* Free surface (or ALE) models. The mesh-coordinate changes now reset the relevant parts of solvers

* New BC interface: add_bc() (#195)

Unifying the two BC functions add_dirichlet and add_natural.
+ The single add_condition() function is much easier to maintain.
+ Improvements in the handling of what 'conds' can be was also made.

 def add_condition(self, f_id, c_type, conds, label, components=None):
        """
        Add a dirichlet or neumann condition to the mesh.

        This function prepares UW data to use PetscDSAddBoundary().

        Parameters
        ----------
        f_id: int
            Index of the solver's field (equation) to apply the condition.
        c_type: string
            BC type. Either dirichlet (essential) or neumann (natural) conditions.
        conds: array_like of floats or a sympy.Matrix
            eg. For a 3D model with an unconstraint x component: (None, 5, 1.2) or sympy.Matrix([sympy.oo, 5, 1.2])
        label: string
            The label name to apply the BC. To find a label/boundary name run something like 
            mesh.view()
        components: array_like, single int value or None.
            (optional) tuple, or int of active conds components to use. Use 'None' for all conds to be used.
            If 'None' and components in 'cond' equal sympy.oo or -sympy.oo those components won't be used.
            eg. For the 3D example cond = (2, 5, 1.2), components = (1,2) the x components is ignored and uncontrainted.
        """

* Cleaning up the warning messages from the tests.

re: no longer needing 'components'

* Adding some info on testing

---------

Co-authored-by: Julian Giordani <julesghub@users.noreply.github.com>
Co-authored-by: Julian Giordani <house.of.jules@gmail.com>

---------

Co-authored-by: Louis Moresi <louis.moresi@anu.edu.au>
Co-authored-by: Ben Knight <55677727+bknight1@users.noreply.github.com>
Co-authored-by: Matthew G. Knepley <knepley@gmail.com>
Co-authored-by: Tyagi <gthyagi@gmail.com>
Co-authored-by: jcgraciosa <12045001+jcgraciosa@users.noreply.github.com>
Co-authored-by: Juan Carlos Graciosa <juan.graciosa@monash.edu>
Co-authored-by: Thyagarajulu Gollapalli <44185316+gthyagi@users.noreply.github.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants