Is it possible to implement a stabilization algorithm for free surface or sticky air models?

[ythaha]

When i test subduction models with low density&viscosity sticky air to simulate a free surface condition, computation seems much more time consuming than without air sometimes, and the velocity and/or temperature in the air show some oscillations around ridge, timestep was limited to much shorter dt( air velocity too large some places).

I know generally when there is an interface with large density contrast, stabilization requires a smaller timestep to avoid sailor drunk effect, but this can be annoying. There is a paper on stabilization algorithm for free surface, where by adding an item on regular FEM formulation the simulation will be more stable —— Kaus, B.J.P., et al., A stabilization algorithm for geodynamic numerical simulations with a free surface. Phys. Earth Planet. In. (2010), doi:10.1016/j.pepi.2010.04.007

So, i'm just curious that whether it is possible to implement this method in Underworld2? In particular, i wonder if it can be achieved by using underworld.systems.sle module? I can avoid some unstable situation in underworld2, so i'm not in an urgent need. I'm just asking to know the possibility and workload for adding this algorithm in uw2. Thanks for any reply!

[julesghub]

Hi @ythaha, We are currently testing a recent implementation of the FSSA algorithm from Kaus et al., the paper you mentioned. Unfortunately the term can't be added by just using the python module underworld.systems.sle - it involves adding some lower level code and extending the python module. This development is taking place on branch jgiordani/fssa-algo.

If you're interested in testing it later I can keep you posted.

[arijitlaik]

@ythaha @julesghub I am interested in this and up for testing, would be great to get this done for the community and "our" models in general

[julesghub]

cheers @arijitlaik - I'll keep that in mind.

[ythaha]

Amazing, that's very helpful, sure i'm interested in this. I'll test it later. Thank you! @julesghub

[rbeucher]

We now have stabilization for sticky air model. We need to add the Free Surface option. We need to extract the normals of the elements. We should be able to leverage PETSc for this. @julesghub is on it.

[LemonBoy68]

Hi~

If I want to view the relevant information of the Le matrix ("lump" matrix, Kaus, 2010), and output it on the screen. I try to add cout (cpp) or PetscPrintf (petsc) and other related command lines in the corresponding cpp file before compiling underworld. After the compilation is successful, I run Benchmark (Kaus, 2010, including free surface) and I don't get any output. . I would be very grateful for any advice.

LYNN

[julesghub]

Hi @LemonBoy68,

We can provide element-wise contributions only as the element matrix is assembled and then applied directly to the global stiffness matrix. So a global "lump" matrix doesn't exist in code. What information to you want to get from the matrix exactly?

On a related note see this repository for further models of the Le stabilisation. https://github.com/underworld-community/topography_drunken_sailor Perhaps @NengLu can also help here.

[LemonBoy68]

@julesghub

Thanks. When will the FSSA2 algorithm be called? I'm interested in how Le fits into the element stiffness matrix, and want to see the value of each calculation, I tried this but it didn't work.

PetscPrintf(PETSC_COMM_WORLD, "fem: %f\n", fem); PetscPrintf(PETSC_COMM_WORLD, "elStiffMat[row][col]: %f\n", elStiffMat[row][col]); ......

In addition, I can't open this URL and can't find relevant information. (https://github.com/underworldcommunity/topography_drunken_sailor)

[NengLu]

Hi @LemonBoy68 The repository is now going public, you can check the implementation of FSSA for the RTI model in Kaus's paper in the repo directory UWG_FSSA/: https://github.com/underworld-community/topography_drunken_sailor/tree/main/UWG_FSSA

It could be called in the UWGeodynamics module directly by: Model.fssa_factor = 1.0

or in the underworld:

fssa_factor = 1.0
fssa = fssa_factor * buoyancyFn * dt * -1.0
stokes_SLE = uw.systems.Stokes(
                velocityField=velocityField,
                pressureField=pressureField,
                conditions=conditions,
                voronoi_swarm=voronoi_swarm,
                _fn_fssa=fssa, 
                fn_viscosity=viscosityFn,
                fn_bodyforce=buoyancyFn)      

And details about the algorithm could be found in function uw.systems.sle.MatrixSurfaceAssemblyTerm_NA__NB__Fn__ni in underworld/systems/_stokes.py:

        if _fn_fssa:
            bgs = uw.swarm.GaussBorderIntegrationSwarm(mesh, particleCount=gaussSwarm._particleCount)
            self._fssa_term = uw.systems.sle.MatrixSurfaceAssemblyTerm_NA__NB__Fn__ni( 
                                            integrationSwarm = bgs, 
                                            assembledObject  = self._kmatrix,
                                            fn               = _fn_fssa,
                                            mesh             = mesh )
        

[LemonBoy68]

@NengLu Thanks a lot, this is very helpful for me. I can now open the link below, which was a weird thing before.