Negative temperature in StickAir Layer.

[tyszwh]

Hi, everyone,

**I run into problems with stickAir layer.

A simple subduction model can run normally with or without stickAir Layer if I use a constant viscosity setting(multicomponet in Material model). But the model is failed to converge when I add a stick_Air Layer above this model if I use viso plastic model in Material model(constant viscosity setting still work).

The main error messeages arised from temperature in stick Air Layer, which decreases until some negative vaules appear. Sometimes the error message is different when you change some boundaries, but it always related to the solve temperature systems.

Could this be a problem with the temperature boundary setting? Why does the temperature of the stick Air Layer could appear negative?

The PRM file Attached below (Corrected File).**

The Error Messages 1

  Postprocessing:
     Writing graphical output: LP_Right_OPEN_Viscous/solution/solution-00000

*** Timestep 1:  t=54031.2 years, dt=54031.2 years
   Solving temperature system... 12 iterations. evrone

--------------------------------------------------------
An error occurred in line <98> of file </home/tysz/888/CIG/aspect/source/material_model/rheology/dislocation_creep.cc> in function
    double aspect::MaterialModel::Rheology::DislocationCreep<dim>::compute_viscosity(double, double, double, unsigned int, const std::vector<double, std::allocator<double> >&, const std::vector<unsigned int>&) const [with int dim = 2]
The violated condition was: 
    viscosity_dislocation > 0.0
Additional information: 
    Negative dislocation viscosity detected. This is unphysical and should not happen. Check for negative parameters. Temperature and pressure are -2.1390754795340299 K, -5125084.3657702813 Pa. 

Stacktrace:
-----------
#0  aspect: aspect::MaterialModel::Rheology::DislocationCreep<2>::compute_viscosity(double, double, double, unsigned int, std::vector<double, std::allocator<double> > const&, std::vector<unsigned int, std::allocator<unsigned int> > const&) const
#1  aspect: aspect::MaterialModel::Rheology::ViscoPlastic<2>::calculate_isostrain_viscosities(aspect::MaterialModel::MaterialModelInputs<2> const&, unsigned int, std::vector<double, std::allocator<double> > const&, std::vector<double, std::allocator<double> > const&, std::vector<unsigned int, std::allocator<unsigned int> > const&) const
#2  aspect: aspect::MaterialModel::ViscoPlastic<2>::evaluate(aspect::MaterialModel::MaterialModelInputs<2> const&, aspect::MaterialModel::MaterialModelOutputs<2>&) const
#3  aspect: aspect::Simulator<2>::local_assemble_advection_system(aspect::Simulator<2>::AdvectionField const&, dealii::Vector<double> const&, dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > const&, aspect::internal::Assembly::Scratch::AdvectionSystem<2>&, aspect::internal::Assembly::CopyData::AdvectionSystem<2>&)
#4  aspect: aspect::Simulator<2>::assemble_advection_system(aspect::Simulator<2>::AdvectionField const&)::{lambda(dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > const&, aspect::internal::Assembly::Scratch::AdvectionSystem<2>&, aspect::internal::Assembly::CopyData::AdvectionSystem<2>&)#1}::operator()(dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > const&, aspect::internal::Assembly::Scratch::AdvectionSystem<2>&, aspect::internal::Assembly::CopyData::AdvectionSystem<2>&) const
#5  aspect: void dealii::WorkStream::run<aspect::Simulator<2>::assemble_advection_system(aspect::Simulator<2>::AdvectionField const&)::{lambda(dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > const&, aspect::internal::Assembly::Scratch::AdvectionSystem<2>&, aspect::internal::Assembly::CopyData::AdvectionSystem<2>&)#1}, aspect::Simulator<2>::assemble_advection_system(aspect::Simulator<2>::AdvectionField const&)::{lambda(aspect::internal::Assembly::CopyData::AdvectionSystem<2> const&)#2}, dealii::FilteredIterator<dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > >, aspect::internal::Assembly::Scratch::AdvectionSystem<2>, aspect::internal::Assembly::CopyData::AdvectionSystem<2> >(dealii::FilteredIterator<dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > > const&, dealii::identity<dealii::FilteredIterator>::type const&, aspect::Simulator<2>::assemble_advection_system(aspect::Simulator<2>::AdvectionField const&)::{lambda(dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false> > const&, aspect::internal::Assembly::Scratch::AdvectionSystem<2>&, aspect::internal::Assembly::CopyData::AdvectionSystem<2>&)#1}, aspect::Simulator<2>::assemble_advection_system(aspect::Simulator<2>::AdvectionField const&)::{lambda(aspect::internal::Assembly::CopyData::AdvectionSystem<2> const&)#2}, aspect::internal::Assembly::Scratch::AdvectionSystem<2> const&, aspect::internal::Assembly::CopyData::AdvectionSystem<2> const&, unsigned int, unsigned int)
#6  aspect: aspect::Simulator<2>::assemble_advection_system(aspect::Simulator<2>::AdvectionField const&)
#7  aspect: aspect::Simulator<2>::assemble_and_solve_composition(bool, std::vector<double, std::allocator<double> >*)
#8  aspect: aspect::Simulator<2>::solve_single_advection_single_stokes()
#9  aspect: aspect::Simulator<2>::solve_timestep()
#10  aspect: aspect::Simulator<2>::run()
#11  aspect: void run_simulator<2>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool)
#12  aspect: main

The Error Message 2

*** Timestep 3:  t=40237.3 years, dt=9431.1 years
   Solving temperature system... retrying linear solve with different preconditioner...
---------------------------------------------------------
TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped
to avoid a possible deadlock.
---------------------------------------------------------

----------------------------------------------------
Exception 'ExcMessage ("The " + solver_name + " did not converge. It reported the following error:\n\n" + exc.what() + "\n The required residual for convergence is: " + std::to_string(solver_controls.front().tolerance()) + ".\n See " + output_filename + " for convergence history.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <387> of file </home/tysz/888/CIG/aspect/source/simulator/solver.cc> in function
    void aspect::{anonymous}::linear_solver_failed(const string&, const string&, const std::vector<dealii::SolverControl>&, const std::exception&)
The violated condition was: 
    false
Additional information: 
    The iterative advection solver did not converge. It reported the following error:

--------------------------------------------------------
An error occurred in line <1091> of file </home/tysz/dealii-candi/deal.II-v9.2.0/include/deal.II/lac/solver_gmres.h> in function
    void dealii::SolverGMRES<VectorType>::solve(const MatrixType&, VectorType&, const VectorType&, const PreconditionerType&) [with MatrixType = dealii::TrilinosWrappers::SparseMatrix; PreconditionerType = dealii::TrilinosWrappers::PreconditionILU; VectorType = dealii::TrilinosWrappers::MPI::Vector]
The violated condition was: 
    iteration_state == SolverControl::success
Additional information: 
Iterative method reported convergence failure in step 1000. The residual in the last step was 1.79764e+16.

This error message can indicate that you have simply not allowed a sufficiently large number of iterations for your iterative solver to converge. This often happens when you increase the size of your problem. In such cases, the last residual will likely still be very small, and you can make the error go away by increasing the allowed number of iterations when setting up the SolverControl object that determines the maximal number of iterations you allow.

The other situation where this error may occur is when your matrix is not invertible (e.g., your matrix has a null-space), or if you try to apply the wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric or not positive definite). In these cases, the residual in the last iteration is likely going to be large.
--------------------------------------------------------

 The required residual for convergence is: 1788759153.808919.
 See RP_Stick_Air_Right_Open_Viscous_2/solver_history.txt for convergence history.
--------------------------------------------------------

Aborting!

simple_2d_cartesian_Stick_Air.txt

[anne-glerum]

Hi @tyszwh, in the attached prm, the Solver scheme is single Advection, single Stokes. Do you always use this scheme with the visco plastic material model? This material model introduces nonlinearities in the rheology that need to be iterated out, using a scheme like single Advection, iterated Stokes for example. Are the velocities at Timestep 0 realistic? Also, a modified thermal conductivity is often used for Sticky Air, so that the temperature in the air remains reasonable.

[tyszwh]

Hi Dr anne @anne-glerum,

Thanks for your suggestions. This test model is also modified from your case. Yes, most of the time I use the single Advection, single Stokes for testing because it is relatively fast. I will try according to your suggestion and report it back.

[tyszwh]

Hi @anne-glerum

I have tested according to your advice. The temperature problem still exists after I modified the thermal conductivity of air to 10 or 20 W/m/k (the other materials are 4 W/m/k) and used an iterated Stokes instead of the single stokes. I also test different velocity boundaries but it seems that the problem is always related to solving the temperature equation.

Is this a matter of resolution?

[anne-glerum]

With the iterated Stokes solver, how many nonlinear iterations were done? Does the velocity field and the viscosity look ok? Are there very high velocities somewhere?

Looking at the prm, the maximum resolution in the model is ~24 km. The sticky-are is a layer of 100 km, where it should be ok, but the sticky-air/crust interface is probably not well resolved. And what about the compositional fields representing the crust? They need a at least a few elements in the vertical direction as well.

I also noticed that set Skip solvers on initial refinement = true. In this case, I'd set it to false, because refinement is dependent on the viscosity, which depends on the solution, so you want to solve.

Then the repetitions in the geometry model do not lead to more or less square elements, but to an aspect ratio of 1:2.6. I'd change it to

    set X extent = 4000e3
    set Y extent = 760e3
    set X repetitions = 5
    set Y repetitions = 1

[tyszwh]

Hi @anne-glerum

I think it is possible, as you said, that the sticky air / crust interface is not well resolved. The C4 reprenst the Ocean Crust whose thickness is about 30km.

Now, I increase the repetitions and refinement value.

The velocity field and viscosity Field at timestep 0 is shown below(single stokes, The velocity boundary did not change) Is it right now? The velocity in airLayer is large especially at right part.

For now, I'm using iterated stokes. But it is really slow, and need more steps to converge(I set the Max nonlinear iterations = 500).

Is it a problem with my velocity boundary setting(same with the prm file)?

It seems that the Relative nonlinear residual (Stokes system) after nonlinear iteration oscillate obvisouly.

It will take some time to iterate out.

[tyszwh]

Hi, @anne-glerum It seems that it is hard to converge (set Nonlinear solver tolerance = 1e-6 and set Max nonlinear iterations = 500). I upload the output file, so you could check the solve history log. LP_RO.zip

[anne-glerum]

Indeed the velocities in the sticky-air are very high. Even though you have an open boundary on the right, the prescribed velocities on the left boundary are far from volume conservative, and that probably makes the problem a lot harder. Can you adjust the prescribed velocities in such a way that the same volume that flows in also flows out? Now there is about 200 km with inflow at 5 cm/yr and 540 km with outflow at 1 cm/yr.

[tyszwh]

Hi, @anne-glerum

I have already adjusted the prescribed velocities yesterday for finding the problem(>560,5cm/yr, <460; -1.087; the stickAir(>660) is 0cm/yr), But it seems that the problem is not in this place (Still unable to converge). In fact, as long as the stick air layer is removed, the model with Nonlinear rheology can run continuously even with stokes. So, perhaps the main problem is about the air layer. Is there any way to set a constant temperature boundary for the entire air layer?

[anne-glerum]

Do you get reasonable velocities in the sticky air layer when the right boundary is tangential slip instead of open?

[tyszwh]

The model can not calculate when the the bottom,top, right boundary is tangential slip. It should be a problem about mass conservation. How can I calculate the inflow and outflow accurately in the prescribed left boundary to ensure mass conservation?

When I set all tangential boundary,The initial velocity field is abnormal and failed to compute with negative dislocation viscosity.

[anne-glerum]

You don't need to worry about mass conservation, just volume conservation, in this incompressible case. So it's like integrating the velocity along the boundary and the total should be 0.

[tyszwh]

Hi, @anne-glerum It seems that there is not a direct method to appoint the temperature of one compositional ( like stcikAir) as a fixed value (or boundary), the temperature in StickAir changes intensly， which causes a 'Negative dislocation viscosity detected' problem of the model (negative parameters such as temperature and pressure) when model include a stikAir layer. Is it right? Do I need to write a plugin?

[anne-glerum]

No, it's not possible to prescribe the temperature for a specific field throughout a model run. @naliboff also had a suggestion, I think it was to specify the Thermal conductivity instead of the Thermal diffusivity in the visco-plastic material model. This is done by setting Define thermal conductivities to true, otherwise the thermal diffusivities are used to compute the thermal conductivity. In the sticky-air, this can lead to issues, probably because of the low density. Do I remember this correctly @naliboff ?

[tyszwh]

Hi, @anne-glerum That is super helpful. It is a problem with the temperature of StickAir, whose value can not be set as a boundary condition during the model run as other codes. And as you said, the model used the thermal diffusivities are used to compute the thermal conductivity even each compositional field was appointed to a fixed value. I have not set 'Define thermal conductivities to true' in the previous test (sorry about that). I think this problem is solved Thanks for your patience @anne-glerum

[anne-glerum]

Glad to hear temperatures behave now! Happy modelling!

[naliboff]

@tyszwh - Sorry to chime in so late and glad to hear using thermal conductivities, rather than thermal diffusivities fixed the issue!

What I've done in both extensional and subduction systems with sticky air layers is to have a high thermal conductivity that supports an initial small temperature gradient in the sticky air layer (e.g., 10 degrees difference between top and bottom). This is probably the most justifiable approach.

A plugin could be developed to set the temperature to a constant value in the sticky air layer, but then one has decide to exactly what constitutes the sticky air boundary. Unless one has an interface tracking scheme that always guarantees a sharp interface, it means you will need to make decisions about what do at gradational interfaces.

A further note of caution about the sticky air layer: if you are using plasticity and examining the evolution of discrete shear bands with Drucker Prager plasticity, all parameters related to the stick-air interface (resolution, material contrast, averaging schemes, etc, etc) will influence the model evolution. I did a lot of testing with stick air for extensional systems and was not comfortable with how much the sticky air layer influenced the fault evolution and broader system evolution. In nearly all cases I now try to use a free surface.

The handling of sticky-air is much easier in codes that are strictly particle-based and also track temperature on the particles. This removes some of the averaging issues you have seen, but noticeably not all of them.

Last, please note that I am not trying to discourage you from using sticky-air. If you are dealing with large-scale subduction that is not highly dependent on the evolution of discrete shear bands, I think it is a reasonable approach. However, I would do very careful testing to see how cell-wise averaging and averaging between fields affects your results.

Perhaps we should a note about all this into the manual?

Cheers, John

[tyszwh]

Hi @naliboff

Thanks for your reply As you said, the handling of the stickAir method is not very easier in current. The jump of temperature, inflow velocity will make the evolution of the model deviate from our expectations. I also use FreeSurface instead of using it for now. However, I will do more tests on the stickAir model(It should be an option when the free surface is hard to work). And I think it will be helpful to add a note if there is a relatively stable method of using the stickAir model.

Here are some other issues that should be addressed that may be helpful. The initial lithostatic pressure boundary is a problem when adding the initial topography in the model (It will help when adding stickAir layer by using initial topography). An option to select a specific flow law for each assigned composition(Is it related to this? @bobmyhill https://github.com/geodynamics/aspect/pull/4239).

In addition, how do melting work together with the viscoplastic material model( I saw a Melt viscoplastic material model pull requested by you several years ago)?