Nan returned from temperature solver with latent heat in heating model

Hello,

We are building some new subduction models in which we want to use the extended Boussinesq Formulation as we have in the past. However, we get an error at the start of the first time-step related to the the advection solver (The initial residual was: nan; see full error below). I’ve found that the issue lies with the “latent heat” option in the Heating Model. If I comment out just this option, the model runs without issue. You can easily replicate the problem by adding the following lines to the convection-box-3d cookbook:

subsection Formulation

set Formulation = custom

set Mass conservation = incompressible

set Temperature equation = reference density profile

set Enable elasticity = false

end

subsection Heating model

set List of model names = adiabatic heating, shear heating, latent heat

end

I also noticed that the heating model subsection of the Parameter Documentation page
no longer lists a “latent heat” option. (only latent heat melt), but the latent_heat.cc and latent_heat.h files are still in the source and include heating_model subdirectory.

Has something changed with way latent heat is treated in the heating model?? I searched on the forum and in the issues on github, but I couldn’t find anything that appeared related to this.

I’m running version 3.1.0-pre with deal.II 9.6.1. The last time I was running models with these same options (and did not have any issues) we were using Aspect version 2.6.0-pre

Thank you for your help,

Magali

Error message:

*** Timestep 0: t=0 seconds, dt=0 seconds

Solving temperature system… retrying linear solve with different preconditioner…

----------------------------------------------------

Exception ‘ExcMessage (exception_message.str())’ on rank 0 on processing:

--------------------------------------------------------

An error occurred in line <2844> of file </quobyte/billengrp/Software/aspect/source/utilities.cc> in function

void aspect::Utilities::throw_linear_solver_failure_exception(const std::string&, const std::string&, const std::vector<dealii::SolverControl>&, const st

d::exception&, MPI_Comm, const std::string&)

Additional information:

The iterative advection solver in Simulator::solve_advection did not

converge.



The initial residual was: nan

The final residual is: nan

The required residual for convergence is: 1.000000e-50

See output-convbox3d-eba/solver_history.txt for the full convergence

history.



The solver reported the following error:



--------------------------------------------------------

An error occurred in line <2054> of file

</quobyte/billengrp/Software/deal.ii/deal.ii-9.6.1-toolchain-gcc-13.2.0-openmpi5.0.5/deal.II-v9.6.1/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 0. The residual

in the last step was nan.

Hi Magali,

I’m using old version, but i guess you need to specify “Melting entropy change“, like this:

subsection Heating model

set List of model names = adiabatic heating, shear heating, latent heat

   subsection Latent heat
       Melting entropy change = - 300 # maybe -500, i don’t know.
  end

end

Because the error NaN and convergence 1.0e-50 indicate something is not initialized or divided by zero.

Hope it helps!

Mingming

Hi Magali,

I was trying to reproduce your issue by adding the lines you posted to the convection-box-3d cookbook. However, I cannot reproduce your error exactly. The model does crash, but only in time step 8, because the temperature becomes basically plus/minus infinity. However, this seems to be related to the adiabatic heating and not the latent heat (i.e., the model behaves the same if I remove latent heat in the heating models, but if I remove adiabatic heating and keep latent heat, it runs past step 8). This also seems logical to me, since convection-box-3d is based on the simple material model, which has no phase transitions, so the latent heat is computed as zero and including it does not make a difference. I tested this with the current development version (3.1.0-pre (main, ac19c42d2) using deal.II 9.7.0).

So if your error still occurs with the current ASPECT version, it’s probably something different than what I was seeing in convection-box-3d, in particular, I assume you are using a different material model. Could you post a minimum example that exactly reproduces the error?

Best wishes,
Juliane

Hi Mingming,

I’m not doing any kind of melting in my models. The latent heat is only related to phase changes.

Thank you for your suggestion.

Magali

Hi Julianne,

I did run this with the convection_box_3d example to test whether this was something specific about our model set-up or something more generic – so that was my minimally simple model that creates the issues. And I get the same error whether I’m using that convection cookbook or something with a different material model and set-up.

I’ve attached the input file and the log files for the convection_box_3d case. I’ve rerun this with debug mode now and it shows a SIGFPE error.

I’ve also attached the input file and log files for another case where we are creating a plume. Neither of these cases has phase transitions in them, so I thought maybe that was an issue, however, my student who is building on Becky Fildes deep earthquake models gets the same error when he tries to rerun one of the models from her paper published last year (those were run using Aspect 2.6). That model has also used the visco plastic material model with elasticity and multiple phase transitions on, and still we get the same error.

It seems like the only thing that is different between our installations is the version of deal.ii. So, I’ll compile version 9.7.1 and see if this makes the issue go away. If this is the case, then I’d be concerned that it is just hiding bug that got exposed with deal.ii version 9.6.1. I’ll let you know that happens.

Magali

(attachments)

eba_convection_box_3d.prm (5.97 KB)
error_convbox3d.log (3.81 KB)
output_convbox3d.log (3.63 KB)
plume2D_simple.prm (6.76 KB)
error_plume2d.log (3.6 KB)
output_plume2d.log (3.6 KB)

Hi Magali,

I ran your eba_convection_box_3d.prm example with the same version of ASPECT I used yesterday and without any problems, i.e. it finishes the one time step the model takes even in debug mode without any errors. From the error message you posted, you’re right, this looks like in the assembly of the temperature equation, something causes a floating point exception when the latent heat heating model is evaluated, and that then causes the nan values if the model is run in release mode (possible some variable is initialized with nan and then not filled with the correct values).

To me, this seems something that is more likely related to ASPECT than deal.ii. So before you reinstall deal.ii: Did you run the latest version of ASPECT, i.e. is your version 3.1.0-pre fairly recent? If not, maybe it would make sense to update that first.

I also looked at the latent heat model and the most recent change I saw came from this PR: add get_required_properties function to heating plugins by jdannberg · Pull Request #6267 · geodynamics/aspect · GitHub (11 months ago)
Basically, in many places, the heating model can tell whoever calls it what kind of inputs are required to compute the heating model outputs, which means that some inputs do not need to be computed (if they’re not needed). I can see how this might cause a problem like yours in case something went wrong (i.e. in case a heating model does not specify all the properties it needs). But that’s all guesswork since for me the model doesn’t produce the error.

Juliane

Hello Juliane,

I compiled the new version of aspect main as mine was from last June, and this indeed fix the problem.

Thanks for your help getting this figured out.

Magali