"geometric average" material averaging crashes in 3D

[jperryhouts]

Maybe just on my machine, but maybe something we should look into before 1.4 release? Can be reproduced with this input:

(works if you change to Dimension=2 or Material averaging = anything else)

 set Dimension                               = 3
 set Output directory                        = build/output

 subsection Geometry model
   set Model name = chunk
   subsection Chunk
     set Chunk inner radius      = 1
     set Chunk outer radius      = 2
     set Chunk minimum latitude  = 37
     set Chunk maximum latitude  = 49
     set Chunk minimum longitude = -132
     set Chunk maximum longitude = -116
   end
 end

 subsection Model settings
   set Fixed temperature boundary indicators     = inner, outer
   set Tangential velocity boundary indicators   = outer
   set Zero velocity boundary indicators         = inner
 end

 subsection Material model
   set Material averaging    = geometric average
   set Model name            = simple
 end

 subsection Initial conditions
   set Model name = adiabatic

   subsection Adiabatic
     set Age bottom boundary layer = 1e6
     set Age top boundary layer = 5e6
     set Amplitude = 100
     set Position = center
     set Radius = 100e3
     set Subadiabaticity = 0
   end
 end

 subsection Boundary temperature model
   set Model name    = initial temperature
 end

 subsection Gravity model
   set Model name    = radial earth-like
 end

 subsection Mesh refinement
   set Initial global refinement             = 3
   set Initial adaptive refinement           = 0
 end

 subsection Termination criteria
   set Termination criteria = end step
   set End step = 0
 end

[bangerth]

On 05/12/2016 01:46 PM, Jonathan Perry-Houts wrote:

Maybe just on my machine, but maybe something we should look into before
1.4 release? Can be reproduced with this input:

Can you also post the error message you get?
Best
W.

---

Wolfgang Bangerth email: bangerth@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[jperryhouts]

Oh, yes, sorry. It's a SolverControl::NoConvergence exception during the temperature advection step, with a NaN residual.

---

Exception on processing:

---

An error occurred in line <446> of file in function
double aspect::Simulator::solve_advection(const aspect::Simulator::AdvectionField&) [with int dim = 3]
The violated condition was:
false
The name and call sequence of the exception was:
ExcMessage (std::string("The iterative advection solver " "did not converge. It reported the following > error:\n\n") + exc.what())
Additional Information:
The iterative advection solver did not converge. It reported the following error:

---

An error occurred in line <1109> of file in function
void dealii::SolverGMRES::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
The name and call sequence of the exception was:
SolverControl::NoConvergence (accumulated_iterations, last_res)
Additional Information:
Iterative method reported convergence failure in step 0. The residual in the last step was -nan.

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.
...

[bangerth]

Hm. I think I've seen something like this as well in a recent case I still have to debug. Does it matter whether you use geometric or harmonic averaging?

[jperryhouts]

Yeah, all other averaging schemes work. It's just geometric averaging for some reason. I poked around in the averaging function, and didn't notice anything obviously different about that case.

[bangerth]

That's strange. If you print out the averaged values on every cell, do you see anything untoward (e.g., NaNs, INFs, etc)?

[tjhei]

I will take a look.

[tjhei]

[0]PETSC ERROR: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero

*** Timestep 0:  t=0 years
   Solving temperature system... 0 iterations.

Program received signal SIGFPE, Arithmetic exception.
aspect::Assemblers::CompleteEquations<3>::local_assemble_stokes_system_incompressible (this=0x154fff0, pressure_scaling=1.4275000000000002e+27, rebuild_stokes_matrix=true, scratch=..., data=...) at /ssd/aspect-local/source/simulator/assembly.cc:1061
1061                        data.local_matrix(i,j) += ( (use_tensor ?
(gdb) p eta
$2 = inf

it turns out the product of the viscosities overflows to infinity because

p values_out
$14 = std::vector of length 27, capacity 27 = {5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24,
 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+2
4, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e
+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.000000000000000
5e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000
005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24, 5.0000000000000005e+24}

:-)