Limit the numeric solutions (v, T) if there are wrong values in some elements

[sibiaoliu]

I met the same issue when I run several different kinds of models: there are some wrong values on a few elements during the calculation, which results in an extremely small timestep or even stop of the simulation. For example, in one of my crustal compression model with a sticky-air, the temperature (or velocity value) can be > 2000 K (or >100 cm/yr) in some elements in the top sticky-air layer, which is clearly artificial value and wrong.

The reason why it has an unnormal value during solving the equations is not clear. I guess that the selection of the iteration method may be a source for causing this problem.

So, I am wondering is there any way that I can add a limit value (e.g., std::numeric_limits::max()?) to the solution after the simulation of one timestep? By doing so, the wrong value can be corrected to speed up the convergence.

Thank you very much for any suggestions in advance!

BTW, how could I label this issue as a question?

[naliboff]

@sibiaoliu thanks for posting the question!

So, I am wondering is there any way that I can add a limit value (e.g., std::numeric_limits::max()?) to the solution after the simulation of one timestep? By doing so, the wrong value can be corrected to speed up the convergence.

In general, I think this approach should be avoided and we should get to the bottom of what is going wrong. I know some groups manually override the temperature in sticky-air layers, but if the sticky-air viscosity is fixed I don't see why this would have an effect on the velocities in that layer.

Can you provide a bit more information about your model setup and how you are handling the sticky-air layer (e.g., all properties)?

This will be a useful discussion for anyone who wants to use the sticky-air approach!

[sibiaoliu]

@naliboff Hi John.

Here I just list one model example that has the issue; there are some other models that do not use sticky air condition and have the same issue. For example, I run one simple chunk model with Gplates surface velocity and it shows the temperature T in some elements is lower than 273K (top boundary temperature, minimum T) or higher than 1837 K (the bottom boundary temperature, maximum T) after a few hundreds of timesteps.

For the sticky-air case, I run a 2d box subduction model with visco-plastic material model. I use the sticky-air condition rather than the free surface due to the mesh distortion in the top boundary element between two plates that stop the model (there is no smoothing function to allow the maxium deformed angle of the surface element yet in Free surface function). This is the same reseaon why we dont use the free surface in the brittle thrust benchmark.

The material properties of the 10-km-thick sticky-air include: the density - 1kg/m3; heat capacity - 1e6 (very higher value than other material); a constant viscosity - 10e19 Pa s; no plasticity (friction angle and cohesion are 0).

In the sticky-air case, here is a snapshot of unnormal RMS velocity in postprocessing after thoundreds of timesteps:

[naliboff]

Hi Sibiao,

Thank you for the detailed reply! Point by point responses below.

Here I just list one model example that has the issue; there are some other models that do not use sticky air condition and have the same issue. For example, I run one simple chunk model with Gplates surface velocity and it shows the temperature T in some elements is lower than 273K (top boundary temperature, minimum T) or higher than 1837 K (the bottom boundary temperature, maximum T) after a few hundreds of timesteps.

This sounds like an overshoot-undershoot issue with the temperature solver. What CFL value are you using? Decreasing the CFL value may help with this issue. Also, are you using the standard temperature solver scheme (e.g., entropy viscosity verse SUPG)?

The material properties of the 10-km-thick sticky-air include: the density - 1kg/m3; heat capacity - 1e6 (very higher value than other material); a constant viscosity - 10e19 Pa s; no plasticity (friction angle and cohesion are 0)

I would actually set the cohesion to something very large like 1e20 Pa to make the sticky air never yields, but I suspect at the selected viscosity (1e19 Pa s) it is not yielding.

My recollection is that when using sticky air in models that also include temperature effects (e.g., not a constant value like in sandbox models) you also need to use a very high thermal conductivity. Are the selected values from a previous study? @anne-glerum @MFraters - Have you done any subduction model tests with sticky-air?

How about we try the sticky-air approach with the continental extension cookbook and see if we can reproduce the observed behavior and/or find a stable combination of parameters for the air layer?

[sibiaoliu]

Hi John,

1. I use the default CFL value as 0.5. Yes, I use the standard temperature solver scheme, i.e., entropy viscosity.

2. The values I use are from my previous models that work well. In the material model of Visco_plastic, we do not set the value of thermal conductivity directly. Instead, we set thermal diffusivity and heat capacity (and reference density) to calculate it. So the thermal conductivity in my model is very high due to a high value of heat capacity I set.

That's a good idea to have a stable sticky-air model of the continental extension/compression cookbook.

Back to my question :) I am interested that where and how I can set a temperature (and velocity) filter to the solution at the end of each time step. Could someone give me some clues? Thanks!

[gassmoeller]

Hello Sibiao, I agree with John that trying to fix the models is almost always the better approach, since applying the filter might just hide the underlying problem. Still if you want to modify the solution after the solver has finished you can do it the following way:

1. The best place might be either in core.cc::run() after the call to the solve_timestep() function, or inside the solve_timestep function once the solver is done.
2. The new solution is stored inside the class member variable solution. It is a block vector that contains all solution fields (velocity, pressure, temperature, all compositions).
3. You can check how such a vector can be modified inside helper_functions.cc::make_pressure_rhs_compatible where we do some modifications to the pressure block (you would want to do something similar to the temperature or velocity block). Start simple with a small modification. You will also want to get a good understanding of what one of those entries means (i.e. what one degree of freedom means, check the introductory text here: https://dealii.org/developer/doxygen/deal.II/step_2.html).

[sibiaoliu]

Hi Rene, Thank you so much for providing the way!

[gassmoeller]

No problem. I am closing the issue for now. Feel free to reopen or open a new one if necessary.