Free-surface problem: cell giving non-positive volume fraction

[kristofporkolab]

Dear All,

I am working on a 2D subduction setup with a continental block in the lower plate, using a free surface, adaptive mesh refinement, and visco-plastic rheology. I have encountered a problem where one of the cells at the top of the model returns a non-positive volume fraction and crashes the simulation:

“An error occurred in line <2727> of file </trinity/opt/apps/software/aspect/aspect-2.0.1/candi/DEAL2/tmp/unpack/deal.II-v9.0.1/source/fe/mapping_q_generic.cc> in function
dealii::CellSimilarity::Similarity dealii::MappingQGeneric<dim, spacedim>::fill_fe_values(const dealii::Triangulation<dim, spacedim>::cell_iterator &, dealii::CellSimilarity::Similarity, const dealii::Quadrature &, const dealii::Mapping<dim, spacedim>::InternalDataBase &, dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim> &) const [with int dim = 2, int spacedim = 2]
The violated condition was:
det > 1e-12*Utilities::fixed_power(cell->diameter()/ std::sqrt(double(dim)))
Additional information:
The image of the mapping applied to cell with center [1.16145e+06 658711] is distorted. The cell geometry or the mapping are invalid, giving a non-positive volume fraction of -42785.9 in quadrature point 3.”

The problematic cell is at or very close to the deforming surface. The attached density plot shows the 1.16e6 x coordinate tick where it occurred.

Below I paste a part of the input file describing global, mesh refinment, and free surface parameters:
set Dimension = 2
set Nonlinear solver scheme = iterated Advection and Newton Stokes
set Nonlinear solver tolerance = 1e-2
set Max nonlinear iterations = 20
set CFL number = 0.05
set Timing output frequency = 20
set Pressure normalization = no

subsection Solver parameters
subsection Stokes solver parameters
set Linear solver tolerance = 1e-3
set Number of cheap Stokes solver steps = 100
end
subsection Newton solver parameters
set Max pre-Newton nonlinear iterations = 1000
end
end

Model geometry (660 km deep aspect ratio 4:1)

subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 4
set Y repetitions = 1
set X extent = 2640e3
set Y extent = 660e3
set X periodic = true
end
end

Mesh refinement specifications:

subsection Mesh refinement
set Initial adaptive refinement = 4
set Initial global refinement = 4
set Time steps between mesh refinement = 5
set Strategy = viscosity, composition, strain rate, temperature
set Refinement fraction = 0.95

subsection Free surface
set Free surface boundary indicators = top
set Free surface stabilization theta = 0.9
set Surface velocity projection = normal
end

The problem occurs where the free surface is rapidly deforming, so perhaps it is related to one of the elements getting “bow-tied” between two re-meshing. Could anyone advise me on how to tackle this issue?

Thank you and best,
Kristof Porkolab

[MFraters]

Hello Kristof,

I also used to have similar problems. If I remember correctly, the solution for me was to switch the Surface velocity projection from normal to vertical. This means that vertices can only more up and down, which should prevent the problem of cells being able to invert.

Best,

Menno

[kristofporkolab]

Hi Menno,

thanks for the reply! I have learnt from the manual that the vertical projection is more stable, therefore that was the first projection I tried. However my continental block is moving long distances laterally, and the vertical projection makes the top of the continent scrape off as it moves (because the fs is projected vertically). So in order to maintain the proper lateral transport of the free surface together with the continent, I had to change it to normal projection. I hoped that the problem can still be solved even with the normal projection.

Cheers, Kristof

[MFraters]

Hello Kristof,

hmm, that is weird. I have never noticed that kind of behavior with the free surface. Material should be able to move along the free surface. Could you show some figures of that so that I can see if I understand correct what you mean? Could you also post the boundary conditions section of the input file which you are using?

Cheers,

Menno

[kristofporkolab]

Hi Menno,

here is the bc section:

Temperature boundary conditions

subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1573
set Top temperature = 273
end
end

Boundary classifications

subsection Boundary velocity model
set Tangential velocity boundary indicators = bottom
end

subsection Nullspace removal
set Remove nullspace = net x translation
end

The uploaded picture shows that in case of vertical projection the free surface stays at the original position of the continent, but since the continent is moving to the right, the upper crust is scraped off. Switching to normal projection solved this problem, but then I got this other bug my original post was about.

Cheers,
Kristof

[MFraters]

What you are fundamentally are trying to do is large deformation together with a normal projection where you move the grid (at least partially) with that deformation. This means that after some time the grid will become very deformed. The only way to deal with this is to re-mesh, which is not something ASPECT is made to do. So I think you should investigate a solution with the vertical projection.

It seems that the material at the surface is somehow not free to move laterally. My suspicion is that the problem is due to the composition boundary conditions which do not move with the continent. It could be something different, but I can’t exclude (or confirm) it from the current current information and figure (would need to see the composition boundary and initial conditions and a bigger figure with grid). But maybe there are some other ideas around what the problem could be.

Cheers,

Menno

[cedrict]

The standard ‘trick’ is to add a little bit of diffusion to the surface, but that is currently not supported by ASPECT. In light thereof, I would for now give up on using a free surface (and resort to free slip), and/or make sure that the subduction channel is well lubricated (low viscous enough) as it looks like the buckling may be caused by a ‘stuck’ subduction.
Sticky air still remains an option.

[jbnaliboff]

Sorry for jumping into the conversation a bit late here. I agree with @cedrict suggestions to make sure the subduction channel is sufficiently lubricated, which may reduce the amount of buckling.

I’m working with another student on a similar problem and we ran into other issues when the initial subduction geometry was not aligned sufficiently to promote subduction. So, another option is to revisit the subduction

A third option is to use the material averaging (e.g., averaging over entire element) option in the material model, which smooths out the viscosity field and reduces some small scale deformation. No idea if this will fix the free surface issue, but it will help improve solver behavior.

The last (and perhaps best) option is to use sticky air, which I am also trying out for subduction problems.

FYI, surface diffusion may be implemented at some point in the next few months, but best to try and other options in the meantime.

[kristofporkolab]

Thank you for the suggestions, I am currently trying to implement them. I can already say that material averaging did result in a slightly smoother geometry, but the error still occurred. I will try the others as well and report back. Thanks again!

[jbnaliboff]

@kristofporkolab - FYI, A group will be working a solution for smoothing out the free surface this week. We will keep you posted as progress is made and time permitting perhaps you can try it out on your model.

[kristofporkolab]

@MFraters - I have been trying to play a bit with the composition boundary conditions, but no success so far. Here is what I have now:

set Number of fields = 3
set Names of fields = oc, uc, lc
end

Spatial domain of different compositional fields

subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function constants = h=660e3, ar=4, dc=4e3, xt=1200e3, wocrhs=300e3,
xcr=800e3, xcl=400e3,
r=6e5, zs=2e5, kappa=1e-6
set Function expression =
if( (h-y)<dc && x>xcr,
1,
if( sqrt((x-xt)(x-xt)+(h-r-y)(h-r-y))<r &&
sqrt((x-xt)(x-xt)+(h-r-y)(h-r-y))>(r-dc) &&
x>xt &&
y>h-zs, 1,
if( x<xcl && (h-y)<dc, 1, 0
)
)
);
if( x>=xcl && x<=xcr && (h-y)<=15e3,
1, 0
);
if( x>=xcl && x<=650e3 && (h-y)<=30e3 && (h-y)>15e3,
1,
if( x>650e3 && x<=xcr && (h-y)>15e3 && y>(h-30e3)+(x-650e3)/10, 1, 0
)
)
end
end

subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end

Temperature boundary conditions

subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1573
set Top temperature = 273
end

Attached a bigger picture with grid. The continental block is initially from x=400e3 to 800e3; as it is transported to the right, becomes more and more tilted/scraped off by the free surface if I have a vertical projection. This does not happen with normal projection.

Cheers, Kristof