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