add examples of usage of GeoParams in basic PT codes

[tduretz]

We can store some example of usage of GeoParams together with simple PT codes.
So far only VE loading and inclusion VE loading work.
The codes are unstable and they need an improved damping strategy. I will take care of that soon.
Would be great to invite @boriskaus to this repository such that he can follow and contribute.
Once we get all basic scripts work (including future VEVP and power-law VE runs) we can think of moving the
scripts to a more appropriate repo (GeoParams?).

[boriskaus]

Sounds good - as to putting this back to the GeoParams repo: I would prefer to have a separate repository and leave GeoParams strictly for the point wise computations. Note that in JustRelax we are collecting pseudotransient solvers (which should work with ParallelStencil and GeoParams). Perhaps an idea to store some solvers there that don’t rely on external packages like PS?

[tduretz]

Great, I've pushed VE benchmark and VE inclusion with improved damping strategy. They should be much better basis for VEVP implementations. They now include design of PT parameters based on Räss et al. (2022). The damping is applied on the residual rather than on flux components.
Let's keep these examples in a separate repository, which we could somehow link to Julia Geodynamics.

[boriskaus]

That should work

[luraess]

Thanks @tduretz for the contribution! I invited @boriskaus to the repo as well.
I agree that we could at some point move these example to a separate repo that may differ from GeoParams (unless we'd like to have a GD example in there). TBD.

[luraess]

@tduretz can you add a .gitignore and include in there .DS_Store and .vscode?

Also, you may want to have to scripts in a scripts folder to keep the src for potentially have the main module if needed at some point in there (e.g. in case of adding unit/reference tests).

[albert-de-montserrat]

This is something that we still need (I think) to add to the rheology kernels, but you should not do v = MatParam[Phasec[i]].CompositeRheology[1] as it is type unstable. This should be replace by phase_viscosity(MatParam, ε̇iic[i], Phasec[i], args) where

@generated function phase_viscosity(v::NTuple{N,Any}, ε̇ii, phase, args) where N
    quote
        Base.@_inline_meta
        Base.@nexprs $N i -> v[i].Phase === phase && return computeViscosity_εII(v[i].CompositeRheology[1], ε̇ii, args)        
    end
end

This results in about x3.5 speedup in the GeoParams version. I will address this issue in GeoParams when I have a spare minute.

[tduretz]

I now get:
ERROR: UndefVarError: DruckerPrager_regularised not defined
Should this be defined somewhere in the script or is it a new GP functionality?

[boriskaus]

that's new in GeoParams; you should link the code with GeoParams#main, or GeoParams#adm-cuda-check if you also want the latest cuda-check turned off. I can also upload my manifest

[luraess]

In the codes, shouldn't one use updated pressure in momentum equation (as "usual")? That would make the scheme slightly more stable.

[tduretz]

Yes

In the codes, shouldn't one use updated pressure in momentum equation (as "usual")? That would make the scheme slightly more stable.

Yes, this can help a little bit. I do this in the latest example with dilation.

[luraess]

What is that status here? Could I merge this upon resolving the conflicts?

[boriskaus]

Hmm, good question. Yes I think it is ok to merge (haven’t run it in a long time, though)

[tduretz]

I have a suggestion. What about we rename the folder examples, which would contain several examples of what one can do with GP. We could have subfolders like '1_flow_laws' (the most basic stuff were discussing about earlier this week, potentially notebooks), '2_deformation_maps' and '3_2D_simulations'. We should still identify which codes in this branch is worth keeping and then move them to '/examples/3_2D_simulations'.

[luraess]

Organising script within distinct examples with respective folders may be a good idea.

[boriskaus]

Good idea!