Free surface velocity BC crashes when large deformation occurs

[QwZhang]

When I run model with large deformation, ASPECT aborts when grid distortion occurs. Here is a simple test case similar to "dam breaking" of viscous materials with free-surface on the top and right walls and the gravity is the only driving force of the flowing, the grids on top right corner distort and ASPECT works only for 35 time steps before it fail to continue.

Another example is that when I run Cedric Thieulot and Anne Glerum's new materials and demos(cookbooks/crust_deformation_2D.prm) with free surface on the top, things go wrong if there are large deformations(e.g., shear dislocation) of geometry and distortion of grids in that area.

I tried strain rate and velocity refinement to overcome this issue, but doesn't help. Did I miss something important?

[QwZhang]

The first run: callape.prm

set Dimension = 2 set Start time = 0 set End time = 1e6 set Use years in output instead of seconds = true set Linear solver tolerance = 1e-6 set Nonlinear solver scheme = iterated Stokes set Nonlinear solver tolerance = 1e-6 set Max nonlinear iterations = 100 set Number of cheap Stokes solver steps = 0 set CFL number = 0.5 set Output directory = output_callapse set Timing output frequency = 1 set Pressure normalization = no

subsection Geometry model set Model name = box subsection Box set X extent = 1e3 set Y extent = 1e3 end end

subsection Model settings set Include adiabatic heating = false set Include shear heating = false set Tangential velocity boundary indicators = left, bottom set Prescribed velocity boundary indicators = set Free surface boundary indicators = top, right end

subsection Material model set Model name = simple subsection Simple model set Reference density = 3300 set Reference specific heat = 1250 set Reference temperature = 0.0 set Thermal conductivity = 1.0 set Thermal expansion coefficient = 4e-5 set Viscosity = 1.e20 end end

subsection Gravity model set Model name = vertical subsection Vertical set Magnitude = 9.81 end end

subsection Boundary temperature model set Model name = box end

subsection Initial conditions set Model name = function subsection Function set Function expression = 0 end end

subsection Mesh refinement set Initial adaptive refinement = 0 set Initial global refinement = 4 set Refinement fraction = 0.95 set Strategy = velocity, strain rate set Coarsening fraction = 0.05 set Time steps between mesh refinement = 0 set Run postprocessors on initial refinement = true end

[ian-r-rose]

Hello Qingwen,

Is the crustal deformation parameter file unchanged from that in the cookbook? That is, are you able to reproduce what is in the manual?

Instabilities like those you have shown are known to arise when the deformation is large, especially near discontinuities or large breaks in slope. You might try playing with the slope refinement indicator. I would like to add some artificial diffusivity to the free surface advection scheme to damp those out when they come up, but haven't really figured out a good way to do that yet.

As for the dam break problem, I fear that it will always be problematic here. The top right cell will essentially always become inverted and bad things will happen. Keep in mind that the free surface plugin was developed for looking at topography over thermal convection problems, and as such, only works for relatively modest deformation.

Best, Ian

On Sun, Jun 7, 2015 at 6:31 AM, Qingwen Zhang notifications@github.com wrote:

The first run: callape.prm

set Dimension = 2 set Start time = 0 set End time = 1e6 set Use years in output instead of seconds = true set Linear solver tolerance = 1e-6 set Nonlinear solver scheme = iterated Stokes set Nonlinear solver tolerance = 1e-6 set Max nonlinear iterations = 100 set Number of cheap Stokes solver steps = 0 set CFL number = 0.5 set Output directory = output_callapse set Timing output frequency = 1 set Pressure normalization = no

subsection Geometry model set Model name = box

subsection Box set X extent = 1e3 set Y extent = 1e3 end end

subsection Model settings set Include adiabatic heating = false set Include shear heating = false set Tangential velocity boundary indicators = left, bottom set Prescribed velocity boundary indicators = set Free surface boundary indicators = top, right end

subsection Material model set Model name = simple subsection Simple model set Reference density = 3300 set Reference specific heat = 1250 set Reference temperature = 0.0 set Thermal conductivity = 1.0 set Thermal expansion coefficient = 4e-5 set Viscosity = 1.e20 end end

subsection Gravity model set Model name = vertical subsection Vertical set Magnitude = 9.81 end end

subsection Boundary temperature model set Model name = box end

subsection Initial conditions set Model name = function

subsection Function set Function expression = 0 end end

subsection Mesh refinement set Initial adaptive refinement = 0 set Initial global refinement = 4 set Refinement fraction = 0.95 set Strategy = velocity, strain rate set Coarsening fraction = 0.05 set Time steps between mesh refinement = 0 set Run postprocessors on initial refinement = true end

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/aspect/issues/526#issuecomment-109756868.

[QwZhang]

Hi Ian, Thank your for your reply. For crustal deformation, I run ASPECT with crust_deformation_2D.prm, the only major difference is I run it with longer End time to let the shear bend dislocates more and then thing go wrong. I added slope refinement as you suggest afterwards but seems helpless:(.

Regards, Qingwen

[bangerth]

@ian-r-rose 's patch in #537 may help with this, but just as a general observation, you will likely always get into trouble if you have large-scale deformation. There will come a point where the topology of the mesh will simply not be able to go along with the deformation any more. For example in the second picture you have above (the Thieulot & Glerum model) you will eventually come to a point where the cells at the shear band will just be too deformed and there is nothing you can do other than doing a complete remeshing -- which the current code cannot do, and I cannot really see how we would support it either. So I encourage you to keep playing with the dynamic geometry feature, but I'll tell you that you will always be able to get to some point where it breaks.

[QwZhang]

@ian-r-rose @ian-r-rose: Thank you for your efforts. I do it when the new commits merge.

[bangerth]

@QwZhang -- have you had an opportunity to try @ian-r-rose 's patch?

[gassmoeller]

It seems we agreed that the free surface implementation is mostly for modest deformation of one surface, and #537 should improve stability for those issues. Since this issue was inactive for a long time I will close it for now. Feel free to reopen when it is still pressing or new information is available.