Dear all,
I am trying to test the cookbook free-surface with the following 'box with lithosphere boundary indicators.
subsection Geometry model
set Model name = box with lithosphere boundary indicators
subsection Box with lithosphere boundary indicators
set X extent = 480.e3
set Y extent = 120.e3
set X repetitions = 4
set Y repetitions = 1
set Y repetitions lithosphere = 1
set Lithospheric thickness = 60e3
end
end
However, the code couldn’t run for more than one timestep, giving out the error of not converging.
Additional information:
The iterative Stokes solver in Simulator::solve_stokes did not
converge.
The initial residual was: 1.869432e+10
The final residual is: 4.168621e+06
The required residual for convergence is: 1.511625e+06
See output1FS+L-V-8-1/solver_history.txt for the full convergence
history.
The solver reported the following error:
--------------------------------------------------------
An error occurred in line <2778> of file
</home/mazq/software/aspect/aspect/source/utilities.cc> in function
void aspect::Utilities::linear_solver_failed(const string&, const
string&, const std::vector<dealii::SolverControl>&, const
std::exception&, ompi_communicator_t* const&, const string&)
The violated condition was:
false
Additional information:
The iterative (top left) solver in BlockSchurPreconditioner::vmult did
not converge.
The required residual for convergence is: 1.027235e-02
The solver reported the following error:
--------------------------------------------------------
An error occurred in line <497> of file
</home/mazq/software/aspect/tmp/unpack/deal.II-v9.4.0/source/lac/trilinos_solver.cc>
in function
void dealii::TrilinosWrappers::SolverBase::do_solve(const
Preconditioner&) [with Preconditioner =
dealii::TrilinosWrappers::PreconditionBase]
The violated condition was:
false
Additional information:
AztecOO::Iterate error code -2: numerical breakdown
Stacktrace:
-----------
I have attached the input parameter file and the log file here. I would really appreciate it if anyone could help me answer whether free surface and ‘box with lithosphere’ boundary indicators could be implemented at the same time without automatically refining the vertical meshing because the vertical resolution became unreasonably high and it takes forever to run a model, which will become a huge problem in 3D models.
log.txt (52.4 KB)
solver_history.txt (3.1 KB)
1FS+L-V_8:1.prm (5.6 KB)
Hi Maziqi,
Thanks for posting this question to the forum. I tried running your attached PRM file with a fairly recent version of main branch in debug, and had a different error occur a bit earlier (see details at the base of the message).
In short, the error is related to significant deformation of elements located at the top of the left/right boundaries in the lithosphere. This can be seen in the attached image.
The root of this issue is that with this geometry model there are now two additional boundaries (left lithosphere, right lithosphere) that are by default set to no-slip boundaries. The following adjustment will set these boundaries to also be free-slip and allows the model to run for millions of years without issue.
subsection Boundary velocity model
set Tangential velocity boundary indicators = left, right, bottom, left lithosphere, right lithosphere
end
For more information on the details of the geometry model, I recommend reviewing the geometry model documentation and the test cases where it is used (see list below):
- tests/lithosphere_boundary_indicator_3d_periodic.prm
- tests/lithosphere_boundary_indicator.prm
- tests/two_merged_boxes_free_surface.prm
- tests/random_perturbation_merged_boxes.prm
- tests/depth_postprocessor_two_merged_boxes.prm
- tests/lithosphere_boundary_indicator_3d.prm
- tests/lithosphere_boundary_indicator_origin.prm
Cheers,
John
*** Timestep 94: t=1251.39 years, dt=13.8296 years
Solving mesh displacement system… 13 iterations.
An error occurred in line <1072> of file </home/naliboff/software/dealii/v9.4.0/install/tmp/unpack/deal.II-v9.4.0/source/fe/mapping_q.cc> in function
dealii::CellSimilarity::Similarity dealii::MappingQ<dim, spacedim>::fill_fe_values(const typename dealii::Triangulation<dim, spacedim>::cell_iterator&, dealii::CellSimilarity::Similarity, const dealii::Quadrature&, const typename dealii::Mapping<dim, spacedim>::InternalDataBase&, dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>&) const [with int dim = 2; int spacedim = 2; typename dealii::Triangulation<dim, spacedim>::cell_iterator = dealii::TriaIterator<dealii::CellAccessor<2, 2> >; typename dealii::Mapping<dim, spacedim>::InternalDataBase = dealii::Mapping<2, 2>::InternalDataBase]
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 [1875 119062] is
distorted. The cell geometry or the mapping are invalid, giving a
non-positive volume fraction of -12015.4 in quadrature point 3.
Hi John,
Thank you very much for your reply. The problem is indeed solved when I set the 4 and 5 boundaries to free slip according to your suggestion. Thanks again!
I do have a follow-up question. When I set up my model like below:
set Model name = box with lithosphere boundary indicators
subsection Box with lithosphere boundary indicators
set Box origin X coordinate = 0
set Box origin Y coordinate = 0
set X extent = 2640e3
set Y extent = 1320e3
set X repetitions = 2
set Y repetitions = 1
set Y repetitions lithosphere = 1
set Lithospheric thickness = 60e3
end
end
The whole box is 2640 * 1320 km, and the lithosphere is 60 km thick, which makes the aspect ratio of the upper lithosphere box 44:1 and the lower 44:21. Due to the fact that the Y repetitions and Y repetitions lithosphere could only be integer above 1, I couldn’t set the upper box resolution to a relatively coarse resolution.
The resolution in the lithosphere is 82.5:3.75 km = 22:1, like this:
However, we don’t want high resolution in the lithosphere because we don’t focus on the lithosphere in our research topic, which will waste our computing time and become a huge problem when in 3D. The reason why we have to add the ‘box with lithosphere’ is that the model couldn’t set up a free surface and the prescribed velocity boundary at the top simultaneously.
Is there any way that I can set up the model with the ‘box with lithosphere’ geometry and free surface without vertically high resolution in the lithosphere? I would greatly appreciate it if you could help me with this.
Cheers
Ziqi
Hi Ziqi,
Glad to hear the proposed solution resolved the issues.
If I understand the question correctly, I think your best option is to do few (or no) global refinements, and then used location based initial adaptive refinement to specify that the mesh should only be refined below a certain depth. The downside to this it that having fewer global refinement levels may affect solver performance.
Cheers,
John
