Negative values of field when modeling subducting slab

[lhy11009]

Hi all

An issue about the DG method on composition fields. I see that the compositional field for the subducting slab in my model has negative values outside of the slab. I have already included:

set Use discontinuous composition discretization = true
subsection Stabilization parameters
    set Use limiter for discontinuous composition solution = true
    set Global composition maximum = 1, 1, 1, 1
    set Global composition minimum = 0, 0, 0, 0
end

But the problem is still there. I have checked the implementation of the 'limiter' for DG method, I think all the values should be reset within the range of [minimum, maximum]. When I plot the figure, I tracked points of negative 'spcrust' value by threshold and plot them with pseudocolor plot of the field itself, the result is: Here the threshold is set to [min, -0.005]. There are points with negative values smaller than -0.2, as shown here:

I want to know if this is something that I should expect from the DG method? Or did I miss something?

[lhy11009]

@mibillen I have opened this issue here, I think you should first take a look. Here I used the default tolerance for solvers. Here is the .prm file:

set Dimension = 2 set Use years in output instead of seconds = true set Start time = 0 set End time = 20e6 set Output directory = output set Additional shared libraries = /home/lochy/software/aspect/plugins/subduction_temperature2d/libsubduction_temperature2d.so, /home/lochy/software/aspect/plugins/prescribe_field/libprescribed_temperature.so set Pressure normalization = surface set Surface pressure = 0 set Adiabatic surface temperature = 273 set Nonlinear solver scheme = single Advection, single Stokes set Prescribe internal temperatures = true

subsection Discretization set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 2 set Use discontinuous composition discretization = true subsection Stabilization parameters set Use limiter for discontinuous composition solution = true set Global composition maximum = 1, 1, 1, 1 set Global composition minimum = 0, 0, 0, 0 end end

subsection Formulation set Formulation = custom set Mass conservation = incompressible set Temperature equation = reference density profile end

subsection Geometry model set Model name = chunk subsection Chunk set Chunk inner radius = 3.481e6 set Chunk outer radius = 6.371e6 set Chunk maximum longitude = 61.0 set Chunk minimum longitude = 0.0 set Longitude repetitions = 2 end end

subsection Mesh refinement set Initial global refinement = 5 set Initial adaptive refinement = 6 set Minimum refinement level = 5 set Strategy = isosurfaces, minimum refinement function, composition approximate gradient set Refinement fraction = 0.7 set Coarsening fraction = 0.2 set Time steps between mesh refinement = 10 set Run postprocessors on initial refinement = true set Skip solvers on initial refinement = true subsection Isosurfaces set Isosurfaces = 11, 12, spcrust: 0.5 | 1.0; 10, 12, spharz: 0.5 | 1.0; 10, 12, opcrust: 0.5 | 1.0; end subsection Minimum refinement function set Coordinate system = spherical set Variable names = r,phi,t set Function constants = Ro=6.371e6, UM=670e3, DD=100e3 set Function expression = ((Ro-r<UM)? \ ((Ro-r<DD)? 8: 6): 0.0) end end

subsection Boundary velocity model set Tangential velocity boundary indicators = west, east, bottom, top end

subsection Initial temperature model set Model name = subduction 2d temperature end

subsection Boundary temperature model set Fixed temperature boundary indicators = bottom, top set List of model names = spherical constant subsection Spherical constant set Inner temperature = 1673 set Outer temperature = 273 end end

subsection Compositional fields set Number of fields = 4 set Names of fields = spcrust,spharz,opcrust,opharz set Compositional field methods = field,field,field,field end

subsection Initial composition model set List of model names = function set List of model operators = add subsection Function set Coordinate system = spherical set Variable names = r, phi set Function constants = PHIC=0.628319, Ro=6.371e6, RC=4.000e+05, \ ST=2.000e+05, DCS=7.500e+03, DHS=3.520e+04, DCO=7.500e+03, DHO=3.520e+04 set Function expression = (((phi<PHIC)&&(Ro-r<=DCS)) ? 1.0 : \ (((phi>=PHIC)&&(Ro-r<=ST)&&(phiRo<=PHICRo+min(RC,sqrt(abs(2.0RC(Ro-r)-(Ro-r)^2.0)))) && \ (RC-sqrt(abs(((phi-PHIC)Ro)^2.0+((Ro-r)-RC)^2.0))<=DCS) ) ? 1 : 0)); \ (((phi<PHIC)&&(Ro-r>DCS)&&(Ro-r<=DHS)) ? 1.0 : \ (((phi>=PHIC)&&(Ro-r<=ST)&&(phiRo<=PHICRo+min(RC,sqrt(abs(2.0RC(Ro-r)-(Ro-r)^2.0)))) \ &&(RC-sqrt(abs(((phi-PHIC)Ro)^2.0+((Ro-r)-RC)^2.0))>DCS) \ &&(RC-sqrt(abs(((phi-PHIC)Ro)^2.0+((Ro-r)-RC)^2.0))<=DHS))? 1.0 : 0));\ (((Ro-r<=DCO)&&(phiRo>PHICRo+min(RC,sqrt(abs(2.0RC(Ro-r)-(Ro-r)^2.0))))) ? 1.0 : 0);\ (((Ro-r>DCO)&&(Ro-r<=DHO)&&(phiRo>PHICRo+min(RC,sqrt(abs(2.0RC*(Ro-r)-(Ro-r)^2.0))))) ? 1.0 : 0) end end

subsection Boundary composition model set Fixed composition boundary indicators = west,east set List of model names = initial composition subsection Initial composition set Minimal composition = 0 set Maximal composition = 1 end end

subsection Material model set Model name = visco plastic set Material averaging = harmonic average subsection Visco Plastic set Reference temperature = 273 set Reference viscosity = 1e20 set Adiabat temperature gradient for viscosity = 7.7075e-09 set Minimum strain rate = 1.e-20 set Reference strain rate = 1.e-15 set Minimum viscosity = 1e19 set Maximum viscosity = 1e24 set Phase transition depths = background:670e3 set Phase transition widths = background:5e3 set Phase transition temperatures = background:1662.0 set Phase transition Clapeyron slopes = background:-2e6 set Thermal diffusivities = 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6,1.0e-6 set Heat capacities = background:1250.0|1250.0, spcrust:1250, spharz:1250, opcrust:1250, opharz:1250 set Densities = background:3300.0|3300.0, spcrust:3000, spharz:3235, opcrust:3000, opharz:3235 set Thermal expansivities = background:3.1e-5|3.1e-5, spcrust:3.1e-5, spharz:3.1e-5, opcrust:3.1e-5, opharz:3.1e-5 set Viscosity averaging scheme = harmonic set Viscous flow law = diffusion set Yield mechanism = drucker set Grain size = 1.0000e-02 set Prefactors for diffusion creep = background:1.4250e-15|1.0657e-18, spcrust:5.0000e-20, spharz:1.4250e-15, opcrust:1.4250e-15, opharz:1.4250e-15 set Grain size exponents for diffusion creep = background:3.0000e+00|3.0000e+00, spcrust:0.0000e+00, spharz:3.0000e+00, opcrust:3.0000e+00, opharz:3.0000e+00 set Activation energies for diffusion creep = background:3.1700e+05|3.1700e+05, spcrust:0.0000e+00, spharz:3.1700e+05, opcrust:3.1700e+05, opharz:3.1700e+05 set Activation volumes for diffusion creep = background:4.0000e-06|1.5000e-06, spcrust:0.0000e+00, spharz:4.0000e-06, opcrust:4.0000e-06, opharz:4.0000e-06 set Prefactors for dislocation creep = background:6.859e-15, spcrust:6.859e-15, spharz:6.859e-15, opcrust:6.859e-15, opharz:6.859e-15 set Stress exponents for dislocation creep = background:3.5, spcrust:3.5, spharz:3.5, opcrust:3.5, opharz:3.5 set Activation energies for dislocation creep = background:480.0e3, spcrust:480.0e3, spharz:480.0e3, opcrust:480.0e3, opharz:480.0e3 set Activation volumes for dislocation creep = background:11.0e-6, spcrust:11.0e-6, spharz:11.0e-6, opcrust:11.0e-6, opharz:11.0e-6 set Angles of internal friction = 25.0, 25.00, 25.00, 25.00, 25.00 set Cohesions = 50.e6, 50.e6, 50.e6, 50.e6, 50.e6 end end

subsection Gravity model set Model name = ascii data end

subsection Prescribed temperatures subsection Indicator function set Coordinate system = spherical set Variable names = r, phi set Function constants = Depth=1.45e5, Width=2.75e5, Ro=6.371e6, PHIM=1.065 set Function expression = (((r>Ro-Depth)&&((rphi<Width)||(r(PHIM-phi)<Width))) ? 1:0) end subsection Temperature function set Coordinate system = spherical set Variable names = r, phi set Function constants = Depth=1.45e5, Width=2.75e5, Ro=6.371e6, PHIM=1.065,\ AGEOP=1.26144e15, TS=2.730e+02, TM=1.673e+03, K=1.000e-06, VSUB=1.58549e-09, PHILIM=1e-6 set Function expression = ((r(PHIM-phi)<Width) ? TS+(TM-TS)(1-erfc(abs(Ro-r)/(2sqrt(KAGEOP)))):\ (phi > PHILIM)? (TS+(TM-TS)(1-erfc(abs(Ro-r)/(2sqrt((KRophi)/VSUB))))): TM) end end

subsection Postprocess set List of postprocessors = visualization, velocity statistics, temperature statistics, depth average subsection Depth average set Number of zones = 50 set Output format = txt set Time between graphical output = 0 end subsection Visualization set List of output variables = density, viscosity, error indicator set Output format = vtu set Time between graphical output = 0.1e6 set Number of grouped files = 0 end end

[mibillen]

@lhy11009 the DG-BP method can only return the exact limits [0 1] if AMR is not used. However, when AMR is used the coarsening operation does not preserve the cell average and this introduces an error, which is itself element size dependent. This is specifically discussed in the paper by He et al., PEPI 2017 (http://dx.doi.org/10.1016/j.pepi.2016.12.001). See section 3.4.2 and 4.0, and Table 2. So we actually expect an overshoot/undershoot error, especially at the boundary of the compositional fields where element size is changing due to AMR. The way to minimize this error is to define the AMR parameters so that the compositional boundaries stay at a refined level, and do not get coarsened and then refined again, repeatedly. The other option is to use a different kind of element (one that is divergence free, such as a Raviart-Thomas element), but this comes at a significant computational cost. So, we need to decide if 1) we can decrease the size of the error with better/different choice of AMR strategies or just a bit more refinement and 2) whether the size of the over/under shoots is important to the problem we are solving. Specifically does it effect the rheology in anyway? Therefore, I do not think there is any issue with the code, but rather we need to improve how we use the AMR and DG-BP methods together.

[lhy11009]

The way to minimize this error is to define the AMR parameters so that the compositional boundaries stay at a refined level, and do not get coarsened and then refined again, repeatedly. Let's try that first. Specifically does it effect the rheology in anyway? I checked and I found that the viscosity structure is not refined well in the computation, let's check that in a later case.

[gassmoeller]

Just adding a third way to solve this with brute-force: You could also limit the input values inside the material model, e.g. instead of in.composition[q][i] you could use std::max(in.composition[q][i], 0.0). This is ok as long as you know that negative values are only produced by numerical over/undershoots and should not influence the solution. Of course you need to make sure that the oscillations do not increase over time (which would no longer be reflected in the material model and could therefore go unnoticed), but that should not be the case here since the dg method should be stable.

[mibillen]

Thanks Rene, This is what we decided to do. -Magali

On Sep 1, 2020, at 1:18 PM, Rene Gassmoeller notifications@github.com wrote:

Just adding a third way to solve this with brute-force: You could also limit the input values inside the material model, e.g. instead of in.composition[q][i] you could use std::max(in.composition[q][i], 0.0). This is ok as long as you know that negative values are only produced by numerical over/undershoots and should not influence the solution. Of course you need to make sure that the oscillations do not increase over time (which would no longer be reflected in the material model and could therefore go unnoticed), but that should not be the case here since the dg method should be stable.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/3822#issuecomment-685108429, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACS64URPLAQB6LK3I524QGTSDVJKBANCNFSM4QIFV7WA.

---

Professor of Geophysics Chair, Geology Graduate Program Geophysics Minor Advisor

Department of Earth and Planetary Sciences Univ. of California, Davis Davis, CA 95616 2129 Earth & Physical Sciences Bldg. Office Phone: (530) 752-4169 https://magalibillen.faculty.ucdavis.edu/

Pronouns: She/her/hers First Generation College Student: https://firstgen.ucdavis.edu/

[maxrudolph]

Hi all - I found this open issue when trying to use the bound preserving limiter in 3D, spherical coordinates on a DG temperature solution. It does not appear to work, and I am seeing min/max temperature values more than 1000K outside of the bounds [0,2500] that I set. I was wondering if the bound preserving limiter has been tested in curvilinear coordinates and whether there is a reason that it might not work?

[naliboff]

@maxrudolph - thanks for bringing this up and unfortunately there was another recent example (see #4131) where the DG limiter was clearly not working in a spherical geometry.

We discussed this in the early summer and the decision was to prevent use of the limiter in spherical geometries until some benchmarks had been done. I simply forgot to make the PR after that meeting, but it is back on the top of my list now.

We discuss reasons why the limiter might need to be modified for spherical geometries, but I don't recall the details. Clearly it is not working in the case presented here and #4131.

Do you see similar temperature overshoots without the limiter?

[maxrudolph]

@naliboff - yes, the issue also occurs without the limiter. I'm not sure that anyone else has tried using DG for temperature in spherical geometry. It produces over/undershoots. I went down this path just trying to find some way to obtain a reasonably accurate solution to the energy equation in spherical geometry. SUPG seems to be plagued by oscillations in spherical geometry as well.

[naliboff]

@maxrudolph - As far as I know, it is only @erinheilman that has tried with the same issues. Do you have a strongly variable viscosity in your simulations. My recollection is that the models presented by @erinheilman showed few to no temperature overshoots as the viscosity contrast (as set by min/max viscosity) decreased. As discussed in #4131, we have some ideas for how to fix this that are currently being tested.

In the meantime, I'm working on a PR to prevent use of the limiters in a spherical geometry (trickier than I thought).

[maxrudolph]

I think that there may be more fundamental issues with the use of DG for the energy equation in spherical coordinates. I saw that Erin's models require an iterative solver. However, all of the models that I've set up have strongly temperature-dependent viscosity but are Newtonian. I have not tried to run an isoviscous case but this wouldn't be too hard to set up. I do not have time between now and July 26 to dive super deep into this, but could do so starting on the 27th.

[naliboff]

Yes, there very well maybe more fundamental issues with DG in spherical geometries. @egpuckett @bangerth @tjhei - Can you comment here?

Max, when you have time I would be really interested to see the following tests:

1. If temperature overshoots of the same magnitude occur without DG
2. How limiting the viscosity affects the temperature overshoots
3. If an iterative advection scheme actually helps here

Independent of your tests, I guess one thing that could be done is redo the spherical benchmarks with DG. @egpuckett - Do recall if any of these benchmarks were performed when DG was implemented?

[maxrudolph]

@naliboff I should mention that the reason that I was experimenting with DG for temperature is that I get overshoots and eventually convergence issues with SUPG. In CitcomS, I think that this is avoided by applying the filter of Lenardic and Kaula (1993). An alternative approach here might be to implement the Lenardic filter in ASPECT. It is at line 689 of this file: https://github.com/geodynamics/citcoms/blob/master/lib/Advection_diffusion.c

[bangerth]

People use DG methods on arbitrary meshes all the time, but one has to be careful in how the penalty parameters are defined in that case because it's no longer so clear what the "mesh width h" should actually be.

My best guess is that this was never tested when @egpuckett and @yinghe616 implemented it back in the day. I think that only @egpuckett would know for sure.

[mibillen]

Hello @maxrudolph, @naliboff, @bangerth , @tjhei, @egpuckett, @lhy11009. I've been looking into the issue the DG method with the limiter because we want to use this for composition (in chunk geometry in 2D and 3D).

I think there's a bit of history that's been lost that might help with tracking down the issue with the DG method. First, I believe that the DG solver that is in Aspect was implemented by Sam Cox (https://github.com/spco), in this pull request: https://github.com/geodynamics/aspect/pull/661 and not by Ying. His focus was on using this for the temperature advection equation and added in the composition because he could.

About the same time is when Ying, Gerry and I were working on this. I think at first we did not know that Sam Cox was also working on this. From the merged pull requests, it seems that Ying only implemented the Bound Preserving limiter in Aspect (https://github.com/geodynamics/aspect/pull/864 and https://github.com/geodynamics/aspect/pull/971) using the DG implementation that was done by Sam Cox.

However, as far as I can tell the bound preserving limiter was only developed and tested for the composition advection equation (without diffusion) in 2d and 3d cartesian. I do not think the limiter as implemented will work on the advection-diffusion equation - there are 2 paragraphs in Ying's paper indicating that this is a future goal and points to 2 recent papers for information on how to do this.

it is still not clear to me whether DG with BP-limiter has any issues with non-cartesian geometries for composition advection. We can put the falling sphere test in a 2D chunk slice and 3D chunk and test if there are crazy overshoots using fixed mesh or refinement only (coarsening will introduce errors because the cell average is not preserved). These tests might help narrow down the source of the errors with the DG method... are they due to geometry or because the DG implementation for the temperature equation has its own problems?