geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
224 stars 235 forks source link

Manifold geometries with a free surface crach #612

Closed MFraters closed 8 years ago

MFraters commented 8 years ago

I am testing the ellipsoidal geometry model, and for this I am testing how it works together with different parts of ASPECT, and especially the free surface. During the testing I stumbled upon a curious bug which seems to be in the constraints part of the free surface code. I have been able to get some constrains on the bug, in the sense that it occurs in combination with at least the spherical shell, chunk and ellipsoidal chunk geometries when using tangential velocity boundary conditions (all the geometries which use manifolds).

The simplest way to reproduce this error is by taking the shell_simple_3d.prm cookbook file and add

set Pressure normalization = no set CFL number = 0.1

and set in the model section:

set Zero velocity boundary indicators = set Tangential velocity boundary indicators = inner,east,south,west set Prescribed velocity boundary indicators = set Free surface boundary indicators = outer

I also had to remove depth average from the list of postprocessors, because it gave some strange error.

The bug shows itself at the beginning of timestep 2:

*** Timestep 2:  t=20425.1 years
------------------------------------------
Got error 2 in row 7537 of proc 15 when trying to add the columns:
481 482 483 484 3052 3053 3054 3055 5194 5195 5196 5197 7537 7538 7539 7540

Matrix row  has the following indices:
7537 7538 7539 7540 7546 7547 7551 7552 7570 7571 7575 7576 7587 7588 7597 7598 7599 7600 7606 7607 7611 7612 7630 7631 7635 7636 7647 7648 482 484 483 480 3053 3051 3055 3054 3101 3103 3102 3099 5233 5194 5196 5245 5197 5244 5232 5228 5195 5227

and in the error output the following:

--------------------------------------------------------
An error occurred in line <1659> of file </home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/source/lac/trilinos_sparse_matrix.cc> in function
    void dealii::TrilinosWrappers::SparseMatrix::add(dealii::TrilinosWrappers::SparseMatrix::size_type, dealii::TrilinosWrappers::SparseMatrix::size_type, const size_type*, const double*, bool, bool)
The violated condition was:
    ierr <= 0
The name and call sequence of the exception was:
    ExcAccessToNonPresentElement(row, col_index_ptr[0])
Additional Information:
You tried to access element (10005/3087) of a sparse matrix, but it appears to not exist in the Trilinos sparsity pattern.

Stacktrace:
-----------
#0  /home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/build/lib/libdeal_II.g.so.8.3.pre: dealii::TrilinosWrappers::SparseMatrix::add(unsigned int, unsigned int, unsigned int const*, double const*, bool, bool)
#1  /home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/build/lib/libdeal_II.g.so.8.3.pre: void dealii::ConstraintMatrix::distribute_local_to_global<dealii::TrilinosWrappers::SparseMatrix, dealii::TrilinosWrappers::MPI::Vector>(dealii::FullMatrix<double> const&, dealii::Vector<double> const&, std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::TrilinosWrappers::SparseMatrix&, dealii::TrilinosWrappers::MPI::Vector&, bool, dealii::internal::bool2type<false>) const
#2  ./aspect: aspect::Simulator<3>::FreeSurfaceHandler::solve_elliptic_problem()
#3  ./aspect: aspect::Simulator<3>::FreeSurfaceHandler::execute()
#4  ./aspect: aspect::Simulator<3>::solve_timestep()
#5  ./aspect: aspect::Simulator<3>::run()
#6  ./aspect: main
--------------------------------------------------------

[node017:08191] *** Process received signal ***
[node017:08191] Signal: Aborted (6)
[node017:08191] Signal code:  (-6)
[node017:08191] [ 0] /lib64/libpthread.so.0[0x3dc6c0f490]
[node017:08191] [ 1] /lib64/libc.so.6(gsignal+0x35)[0x3dc6832945]
[node017:08191] [ 2] /lib64/libc.so.6(abort+0x175)[0x3dc6834125]
[node017:08191] [ 3] /home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/build/lib/libdeal_II.g.so.8.3.pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseEb+0x3b)[0x7ffff5dc5f38]
[node017:08191] [ 4] /home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/build/lib/libdeal_II.g.so.8.3.pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_16TrilinosWrappers12SparseMatrix28ExcAccessToNonPresentElementEEEvNS1_17ExceptionHandlingEPKciS8_S8_S8_T_+0x34)[0x7ffff5c770e6]
[node017:08191] [ 5] /home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/build/lib/libdeal_II.g.so.8.3.pre(_ZN6dealii16TrilinosWrappers12SparseMatrix3addEjjPKjPKdbb+0xac7)[0x7ffff5c73231]
[node017:08191] [ 6] /home/aspect/deal.ii/deal.ii_2015-06-04-v11.4/dealii/build/lib/libdeal_II.g.so.8.3.pre(_ZNK6dealii16ConstraintMatrix26distribute_local_to_globalINS_16TrilinosWrappers12SparseMatrixENS2_3MPI6VectorEEEvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKSt6vectorIjSaIjEERT_RT0_bNS_8internal9bool2typeILb0EEE+0x5bf)[0x7ffff5783da3]
[node017:08191] [ 7] ./aspect(_ZN6aspect9SimulatorILi3EE18FreeSurfaceHandler22solve_elliptic_problemEv+0x6de)[0x906d08]
[node017:08191] [ 8] ./aspect(_ZN6aspect9SimulatorILi3EE18FreeSurfaceHandler7executeEv+0x61)[0x9083ab]
[node017:08191] [ 9] ./aspect(_ZN6aspect9SimulatorILi3EE14solve_timestepEv+0x11a)[0x8eb0f2]
[node017:08191] [10] ./aspect(_ZN6aspect9SimulatorILi3EE3runEv+0x44f)[0x8f3d9b]
[node017:08191] [11] ./aspect(main+0x4c5)[0x99a751]
[node017:08191] [12] /lib64/libc.so.6(__libc_start_main+0xfd)[0x3dc681ec9d]
[node017:08191] [13] ./aspect[0x625c89]

I am not exactly causes the problem. Could it for example be caused by the decoupling of the manifold by the free surface code? If needed I can also attach input files where the problem occurs for both the chunk and ellipsoidal chunk geometry.

Many thanks!

Menno

ian-r-rose commented 8 years ago

Hi Menno,

This does indeed seem like a bug. I am having a bit of a hard time reproducing the error here -- can you confirm that it only happens when there are free-slip boundary conditions? Also, any chance for a prm file that I can run on my workstation (preferably in 2D)?

Ian

MFraters commented 8 years ago

Hey Ian, Thanks for your quick reply. The 2d case is a lot faster (cookbooks/shell_simple_2d.prm), and also shows the bug. Just change the same things as above (but you need left and right in stead of east,west and south). I forgot to say that in the 3d case you also need to add set Opening angle = 90. In the 2d case I can say that the problem only occurs when the inner boundary is set to tangential, left and right can be either zero or tangential without a problem.

ian-r-rose commented 8 years ago

Hmm, I just tried shell_simple_2d.prm, and it seems to work okay on my machine. Care to attach your prm with the appropriate modifications so we can be 100% sure that we are doing the same thing?

I am on the ASPECT master branch, using deal.ii 8.3.

MFraters commented 8 years ago

I have noticed a difference in behaviour between my laptop (with an oder version of deal.ii) and the cluster with a recent version (june 2015 from github). My laptop does generate the error with the 2d version, but the cluster doesn't. But the 3d version crashed in both cases. The problem is dependent on the resolution which is used, so this is the smallest 3d case I found (of I take a lower adaptive refinement the problem disappears):

# A simple setup for convection in a 3d shell. See the
# manual for more information.

set Dimension                              = 3
set Use years in output instead of seconds = true
set End time                               = 1.5e9
set Output directory                       = output
set Pressure normalization = no
set CFL number = 0.1

subsection Material model
  set Model name = simple

  subsection Simple model
    set Thermal expansion coefficient = 4e-5
    set Viscosity                     = 1e22
  end
end

subsection Geometry model
  set Model name = spherical shell

  subsection Spherical shell
    set Opening angle = 90
    set Inner radius  = 3481000
    set Outer radius  = 6336000
  end
end

subsection Model settings
  set Zero velocity boundary indicators       =
  set Tangential velocity boundary indicators = inner,east,west,south
  set Prescribed velocity boundary indicators =

  set Fixed temperature boundary indicators   = inner, outer
 set Free surface boundary indicators = outer
  set Include shear heating                   = false
end

subsection Boundary temperature model
  set Model name = spherical constant
  subsection Spherical constant
    set Inner temperature = 1973
    set Outer temperature = 973
  end
end

subsection Initial conditions
  set Model name = function
  subsection Function
    set Function expression = 1.473e3
  end
end

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

subsection Mesh refinement
  set Initial global refinement          = 2
  set Initial adaptive refinement        = 2
  set Strategy                           = temperature
  set Time steps between mesh refinement = 15
end

subsection Postprocess
  set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics #, depth average

  subsection Visualization
    set Output format                 = vtu
    set Time between graphical output = 1e6
    set Number of grouped files       = 1
  end

  subsection Depth average
    set Time between graphical output = 1.5e6
    set Output format                 = vtu
  end
end

subsection Checkpointing
  set Steps between checkpoint = 50
end