Temperature solution issues with 3D spherical geometry

[erinheilman]

I've been seeing issues while running 3d spherical shell with 90 degree opening models. In a purely isoviscous case where the boundary temperatures are prescribed as 273 on the top and 2000 on the bottom, the temperature solution runs out of control to numbers way out of bounds. I'm not sure if it is a solver issue or a geometry issue, this model runs in debug mode with no errors as well. I have also discussed this with @gassmoeller previously. 3D_sphere90_01.txt

log_3D_sphere_90_01.txt

[bobmyhill]

A ?related problem... when I change the Boundary temperature model in @erinheilman's input file to "Constant" and assign values to the top and bottom boundaries like so:

subsection Boundary temperature model set List of model names = constant subsection Constant set Boundary indicator to temperature mappings = bottom:2000, top:273 end end

ASPECT doesn't fix the temperature at those values. This isn't limited to the 90 degree opening; the same thing happens for the complete spherical shell.

[Follow up: Actually, this is an entirely unrelated issue which results from not setting "Fixed temperature boundary indicators". I'll open a new PR fixing this]

[erinheilman]

I've also had this happen in a 3D chunk geometry model, so it may just be a 3D spherical issue

[anne-glerum]

Does it still happen when you prescribe boundary velocity conditions on the radial boundaries?

[bobmyhill]

Yes, it does. I've been playing with this for a few hours, and am stumped.

[tjhei]

I have an idea. It looks like the boundary ids are not assigned correctly. I will check.

[tjhei]

I can not reproduce this on my machine with deal.II 9.2. Trying 9.3 next...

[erinheilman]

I have deal.II 9.3.0 on my computer when I ran the model, but I had still had problems with this when using deal.II 9.2

[tjhei]

Do you also see the problem with "global refinement = 2" ?

@zjiaqi2018 can you try to reproduce the problem on your system?

[bobmyhill]

I see the problem with global refinement = 1 or 2. It looks like a quadratic over/undershoot to me, but why it should do that is beyond me.

[erinheilman]

I've tried running models with higher mesh refinement in both AMR and finite mesh and that will reduce the magnitude of the over/undershoots but the temperature will still go out of bounds

[tjhei]

Is this visible in the first output at time 0?

[erinheilman]

Typically no, it will take a few timesteps for the temperature to go out of bounds

[bobmyhill]

At global refinement = 1, it takes ~20 time steps to reach steady state.

[tjhei]

ohhh, ok

[bangerth]

Is this a case where the model just develops a temperature boundary layer and the mesh is simply too coarse to resolve it and then you end up with numerical oscillations that manifest as over/undershoots? It would be interesting to visualize on a finer mesh than the computational mesh, just to get a better idea of how the solution really looks. Are you setting "Interpolate output" to "true" in your .prm file?

[tjhei]

Yes, I can see over/undershoots with refine=2, but they are not on the boundary itself but one cell inside.

[tjhei]

Erin, do you agree that this is just an issue with resolution?

[erinheilman]

I agree that the problem is exacerbated by having a lower resolution mesh, but I still get temperature over/undershoots when using higher resolution mesh (both with finite and amr). For example, I ran this 2D cylindrical model this week. The temps at the boundaries are set to 273 on the outer and 2973 on the inner and it has high refinement at the boundaries but the temperature still exceeds the boundaries by a lot.

[bangerth]

Can you calculate what the expected width of the boundary layer should be for the values of the thermal conductivity, and how that compares to the mesh size? That is, in essence, what determines whether or not you have over/undershoots: whether you can resolve the boundary layer or not.

[erinheilman]

The expected width of the top thermal boundary layer in this model is ~79 km, and there are 4 mesh cells that span that distance. For reference,

[naliboff]

@erinheilman - can add a screenshot for this model where the temperature is exceeding the expected value? Here it appears the values are roughly correct?

[erinheilman]

Sure, this is an earlier timestep of the same model where the temperatures exceed the boundaries, and the corresponding mesh.

[naliboff]

Summary of a conversation I just had with @erinheilman:

• I think part of the issue is that these models require an iterative scheme for both the Stokes and advection systems.
• My guess is that overshoots in temperature lead to unstable/unrealistic viscosities, which then feed back into the temperature through advection.
• @erinheilman is going to send me a parameter file and I will make some suggestions here (@tjhei - this may be duplicative with what you are already working on)?

[tjhei]

• (@tjhei - this may be duplicative with what you are already working on)?

I am not working on anything related to this PR. Let me know if you want me to look at something...

[naliboff]

@erinheilman - thanks for passing on the parameter file! Below are my suggestions for changing a wide range of parameters.

I think the place to start is with a fairly stable model (low viscosity contrasts, small time steps, conservative solver tolerances) to see if that produces realistic temperature values. If this is the case, then the next step would be to make the problem more realistic (e.g., harder to solve) piece by piece.

As mentioned, you would need to remove the strain weakening/healing from your model to make the iterated Advection schemes work correctly. Hopefully there will be a fix for this soon, but you can always use a "single Advection" scheme if that does not produce the same issues you saw above (but again, start with an iterated scheme).

If all goes well with the suggestions below, it would be worth discussing with @MFraters how to use the Newton solver for your models. Given you are not using pressure dependent plasticity, in theory the Newton solver should be well suited for this problem.

Please let us know how it goes with the suggestions below!

Nonlinear solver scheme and tolerance:

# "single Advection, iterated defect correction Stokes" will work with strain healing/weakening
set Nonlinear solver scheme                = iterated Advection and defect correction Stokes
set Nonlinear solver tolerance             = 1e-5 # you could go down to 1e-4 as well

Time stepping:

# decreasing this should help with advection overshoots, can try increasing back to 0.5 later
set CFL number                             = 0.2 

Solvers

subsection Solver parameters
  subsection Stokes solver parameters
    # The linear solver tolerance below is only used once with Eisentstat walker iterations
    set Linear solver tolerance                        = 1.0e-7
    set Stokes solver type                                     = block GMG  
    set Number of cheap Stokes solver steps     = 200
  end
  subsection Newton solver parameters
     set Maximum linear Stokes solver tolerance   = 1e-2
     set Use Eisenstat Walker method for Picard iterations = true
  end
  #set Temperature solver tolerance        = 1e-4 # delete this line (e.g., use default value)
  #set Composition solver tolerance        = 1e-4 # delete this line (e.g., use default value)
end

Material model (only listing some parameters)

subsection Material model

 # This parameter determines how values are averaged over the entire cell.
 # This will reduce accuracy, but also produce more stable convergence behavior.
 # harmonic average will provide the most "smoothing" 
 set Material averaging = harmonic average

  set Model name = visco plastic
  subsection Visco Plastic

    # reference stuff
    set Reference temperature = 1523

    set Reference viscosity   = 1e21

    set Minimum strain rate   = 1.e-20

    # After ensuring the temperature and viscosities are reasonable, start increasing the 
    # viscosity contrast through the min/max values.
    set Minimum viscosity     = 1e20 # start with high value to produce stable convergence
    set Maximum viscosity     = 1e24 # start with low value to produce stable convergne
end