Still an issue with writing out large files?

[elodie-kendall]

Hi,

I am running some big global 3D shell models (which ran fine in 2D) which are running fine until a certain time step (966) and then the code resumes from Timestep 900 snapshot. This then happens continually and the model can't progress. All parameters look fine so maybe it is that the output files are becoming to large again. I am using the latest deal version (10.0-pre) with the latest ASPECT branch which should have fixed the by the >40GB issue raised earlier. The largest file in my output folder is restart.mesh_fixed.data, at 43GB. I attach my prm and log file.

Thanks, Elodie log.txt

#!/usr/bin/env python3

####  Global parameters
set Dimension                              = 3
set Use years in output instead of seconds = true
set Start time                             = 0
set End time                               = 300.e6  
set Output directory                       = output
set Nonlinear solver scheme                = single Advection, iterated defect correction Stokes
set Nonlinear solver tolerance             = 5e-4
set Max nonlinear iterations               = 150
set Max nonlinear iterations in pre-refinement = 0
set Pressure normalization                 = surface
set Maximum time step                      = 1e6
#set Number of cheap Stokes solver steps    = 200  # aspect 2.2
set Adiabatic surface temperature          = 1700 
set CFL number                             = 1

set Resume computation          = auto 
subsection Checkpointing
  set Steps between checkpoint  = 100
  set Time between checkpoint   = 0
end
subsection Termination criteria
  set Checkpoint on termination           = true
  set End step                            = 1
  set Termination criteria                = end time #end step
end

subsection Solver parameters
  subsection Stokes solver parameters
    set GMRES solver restart length = 250
    set Linear solver tolerance             = 1e-2
    set Maximum number of expensive Stokes solver steps = 1000
    set Number of cheap Stokes solver steps    = 200 
  end
end

#### Model Geometry
subsection Geometry model
  set Model name = spherical shell

  subsection Spherical shell
    set Inner radius  = 3485e3
    set Outer radius  = 6371e3 #5961e3
    set Opening angle = 360
  end
end

#### Mesh refinement
subsection Mesh refinement
#  set Coarsening fraction = 0.2
#  set Refinement fraction = 0.8
  set Initial adaptive refinement        = 3
  set Initial global refinement          = 4
  set Time steps between mesh refinement = 0
  set Strategy                           = minimum refinement function
  subsection Minimum refinement function
    set Variable names      = x,y,z
    set Function constants  = R=6371.e3
    set Coordinate system = cartesian
    set Function expression = if(sqrt(x*x+y*y+z*z)>=R-200e3, 7,if(sqrt(x*x+y*y+z*z)>=R-620e3, 5,if(sqrt(x*x+y*y+z*z)>=R-680e3,5,if(sqrt(x*x+y*y+z*z)<=R-2586e3,5,4)))) 
    #if(depth<560e3,0,if(depth>760e3,0,6))
#        if(depth<560e3,5,if(depth>760e3,5,6))
  end
end

#### Set velocity model, zero bottom, tangential others
#subsection Boundary velocity model
#  set Zero velocity boundary indicators       = inner
#  set Tangential velocity boundary indicators = left, right
#end

#subsection Model settings
#  set Fixed composition boundary indicators   = bottom, top 
#  set Fixed temperature boundary indicators   = top, bottom
#  set Zero velocity boundary indicators       = inner
#  set Tangential velocity boundary indicators = left, right
#  #set Prescribed velocity boundary indicators = 3: function
#  set Free surface boundary indicators = top
#end

#subsection Boundary composition model
  #set Allow fixed composition on outflow boundaries = false  # false for models without melt
#  set Fixed composition boundary indicators         = bottom, top
#end
subsection Boundary temperature model
  #set Allow fixed temperature on outflow boundaries = true  # default
  set Fixed temperature boundary indicators         = bottom, top
end
subsection Boundary velocity model
  set Tangential velocity boundary indicators       = inner, outer
 # set Zero velocity boundary indicators             = outer
end   

# Free surface, where exaggeration -> vertical
# subsection Mesh deformation
#   set Mesh deformation boundary indicators        = top: free surface
  #set Free surface boundary indicators = top
#   subsection Free surface
#     set Surface velocity projection             = vertical
#     set Free surface stabilization theta        = 0.5
#   end
# end

subsection Nullspace removal
   set Remove nullspace = net rotation
end

#subsection Mesh deformation
#  set Mesh deformation boundary indicators = top: free surface
#  subsection Free surface
#    set Free surface stabilization theta        = 0.5
#    set Surface velocity projection             = vertical
#  end
#end

# Number and names of compositional fields
#subsection Compositional fields
#  set Number of fields = 1
#  set Names of fields = plastic_strain #add strain field
#end

# Spatial domain of different compositional fields
# Horizontal layers
#subsection Initial composition model
#  set Model name = function
#  subsection Function
#    set Coordinate system    = depth
#    set Variable names       = depth,t                
#    set Function expression  = \
#                               if(depth>=660e3, 1, 0); \
#                               if(depth<660e3, 1, 0)

#subsection Initial composition model
#  set Model name = function
#  subsection Function
#    set Coordinate system    = cartesian
#    set Variable names       = x,y,z 
#    set Function constants   = R=6371.e3               
#    set Function expression  = \
#                               if(sqrt(pow(x,2)+pow(y,2))<660e3,1,0); \
#                               \
#                               if(sqrt(pow(x,2)+pow(y,2))>=660e3,1,0)#,\ 
##                               if((sqrt(pow(x,2)+pow(y,2))>R-200.e3 && \
##                                  sqrt(pow(x,2)+pow(y,2))<R- 20.e3  && \
##                                  atan(y/x)> 139*pi/180 && \
##                                  atan(y/x)< 141*pi/180) || \
##                                  (sqrt(pow(x,2)+pow(y,2))>R-200.e3 && \
##                                  sqrt(pow(x,2)+pow(y,2))<R- 20.e3 && \
##                                  atan(y/x)> 45*pi/180  && \
##                                  atan(y/x)< 46*pi/180)  || \
##                                  (sqrt(pow(x,2)+pow(y,2))>R-200.e3 && \
##                                  sqrt(pow(x,2)+pow(y,2))<R- 20.e3 && \
##                                  atan(y/x)> 88*pi/180  && \
##                                  atan(y/x)< 92*pi/180)    \
##                                  ,0,1)
#                              #if(sqrt(pow(x-4.e6,2)+pow(y-4.e6,2))<600.e3,0,1)#, \
#                                
##                               if(depth<760e3 && depth>=560e3,-(depth-560e3)/(560e3-760e3),if(depth>760e3, 1, 0)); \
##                               if(depth>560e3 && depth<=760e3,-(depth-760e3)/(760e3-560e3),if(depth<560e3, 1, 0))
#                               # Transitional 
##                               if(depth<760e3 && depth>=560e3,-(depth-560e3)/(560e3-760e3),if(depth>760e3, 1, 0)); \
##                               if(depth>560e3 && depth<=760e3,-(depth-760e3)/(760e3-560e3),if(depth<560e3, 1, 0))
#  end
#end

# Composition fixed on bottom, free on sides and top
#subsection Boundary composition model
#  #set Fixed composition boundary indicators = bottom, top 
#  set List of model names = computeprofile # initial compositionS
#  set Model name = initial composition
#end

##### Set temperature Model
subsection Boundary temperature model
  set Fixed temperature boundary indicators = top, bottom
  set List of model names = spherical constant

  subsection Spherical constant
    set Inner temperature = 4000
    set Outer temperature = 273 #273 for surface ############
  end
end    

#### Temperatur model
subsection Initial temperature model
  #set Model name = adiabatic
  set List of model names =  adiabatic, function
  set List of model operators = add
#  set Model name = adiabatic
  subsection Adiabatic
    set Age bottom boundary layer = 100.e6
    set Age top boundary layer = 100.e6
    set Radius = 0
    set Amplitude = 0
    set Position = center
    set Subadiabaticity = 0
#    subsection Function         
#      set Variable names      = x,y
#      set Function constants  = R=6371.e3
#      set Function expression = \
#                                if(sqrt(x*x+y*y)<=R-660.e3, 1, 0); \
#                                if(sqrt(x*x+y*y)>R-660.e3, 1, 0)
#   end
 end 
 subsection Function
     set Coordinate system = cartesian
     set Variable names      = x,y,z
     set Function constants  = R=6371.e3
     set Function expression = \
                               if(sqrt(x*x+y*y+z*z)<=R, rand_seed(1) * 100 -50 , 0)
 end 

#  subsection Function
#      set Coordinate system = cartesian
#      set Variable names = x,y,z
#      set Function expression = 1
#  end

#  subsection Function
  #  set Coordinate system = depth
  #  set Variable names = depth, t
  #  set Function expression = 1600+0.25*depth/1000
#  end
end

subsection Adiabatic conditions model  # need to initialize for initial adiabatic temp.
  #set Model name = compute profile
  #set Model name = initial profile
  #subsection Compute profile #Initial profile
  #    set Variable names      = x,y
  #    set Function constants  = R=6371.e3
  #    set Function expression = \
  #                              if(sqrt(x*x+y*y)<=R-660.e3, 1, 0); \
  #                             if(sqrt(x*x+y*y)>R-660.e3, 1, 0)
#                                if(sqrt(x*x+y*y)>=R-660.e3, 1, 0); \
#                                if(sqrt(x*x+y*y)<R-660.e3, 1, 0)
  #send
  set Model name = ascii data
  subsection Ascii data model
    set Data directory = $ASPECT_SOURCE_DIR/data/adiabatic-conditions/ascii-data/
    set Data file name = example_isentrope_1700.txt
  end
end

subsection Heating model
  set List of model names = adiabatic heating, shear heating, radioactive decay
  subsection  Adiabatic heating
    set Use simplified adiabatic heating = true
  end
  subsection Radioactive decay
    set Crust defined by composition = false 
    set Crust depth = 10e3
    set Half decay times = 4.468e9, 0.704e9, 14e9, 1.248e9
    set Heating rates = 9.4946e-5, 5.68e-4, 2.6368e-5, 2.8761e-5
    set Initial concentrations crust = 0.0, 0.0, 0.0, 0.0
    set Initial concentrations mantle =0.0082, 0.000059, 0.024, 0.013
    set Number of elements = 4
  end
  #set List of model names = compositional heating
end

############################ MATERIAL MODEL ######################################################################

subsection Material model
  set Model name = compositing #visco plastic #melt simple
############################ MATERIAL MODEL ######################################################################
  subsection Compositing
    set Viscosity = visco plastic
    set Density = ascii reference profile
    set Thermal expansion coefficient = ascii reference profile
    set Specific heat = ascii reference profile
    set Thermal conductivity = ascii reference profile
    set Compressibility = ascii reference profile
    set Entropy derivative pressure = ascii reference profile
    set Entropy derivative temperature = ascii reference profile
    set Reaction terms = visco plastic
  end
    subsection Ascii reference profile
        set Thermal conductivity =  3.0
        subsection Ascii data model
          set Data file name = example_isentrope_1700.txt
          set Data directory = $ASPECT_SOURCE_DIR/data/adiabatic-conditions/ascii-data/
      end
    end

    subsection Visco Plastic 

      set Define transition by depth instead of pressure  = true
      set Phase transition Clapeyron slopes = background: 3.0e6|-2.5e6|9.e6 
      set Phase transition depths           = background: 410e3|660e3|2740e3
      set Phase transition widths           = background: 12e3|23e3|11e3
      set Phase transition temperatures     = background: 1770|1852|2591

      # Reference temperature and viscosity
      set Reference temperature = 293
      set Reference viscosity = 6e22

      # The minimum strain-rate helps limit large viscosities values that arise
      # as the strain-rate approaches zero.
      # The reference strain-rate is used on the first non-linear iteration
      # of the first time step when the velocity has not been determined yet. 
      set Minimum strain rate = 1.e-20
      set Reference strain rate = 1.e-15

      # Limit the viscosity with minimum and maximum values
      set Minimum viscosity = 1e21
      set Maximum viscosity = 1e24

      # Harmonic viscosity averaging
      set Viscosity averaging scheme = harmonic

      # Choose to have the viscosity (pre-yield) follow a dislocation
      # diffusion or composite flow law.  Here, dislocation is selected
      # so no need to specify diffusion creep parameters below, which are
      # only used if "diffusion" or "composite" option is selected.
      set Viscous flow law = composite

#  from cizkov
#     order:                           background    ,        mantle_low        ,        mantle_up
#    set Thermal diffusivities =   8.384146341463415e-07,     8.384146341463415e-07   , 7.309941520467837e-07

    set Thermal diffusivities = 1.e-6
    set Thermal expansivities = background: 3.e-5|2.e-5|2.e-5|2.e-5#, plastic_strain: 3.e-5
    set Heat capacities       =        1250     # specific heat
    set Densities             = background: 3416|3800|4750|4750#, plastic_strain: 3416
#    set Thermal expansivities =      3.0e-5  ,              3.0e-5,                     2.7e-5    
#     Dislocation creep                                                 #(1/2)^3.5 * 1e-17 
    set Prefactors for dislocation creep          = background: 4.54e-15|4.54e-15|1e-30|1e-30#, plastic_strain:1e-22
    set Stress exponents for dislocation creep    = 1.0
    set Activation energies for dislocation creep = background: 1.51e5|1.51e5|4.8e5|4.8e5#, plastic_strain:4.8e5
    set Activation volumes for dislocation creep  = background: 17.e-6|17.e-6|11.e-6|11.e-6#, plastic_strain: 11.e-6
#     Diffusion creep
    set Prefactors for diffusion creep            = background: 1.e-9|1.e-9|1.0e-16|1.0e-16#, plastic_strain: 1.e-9
    set Stress exponents for diffusion creep      = 1.0                                                  #### doesnt exist anymore??? ####
    set Activation energies for diffusion creep   = background: 3.35e5|3.35e5|2.e5|2.e5#, plastic_strain: 3.35e5
    set Activation volumes for diffusion creep    = background: 4.8e-6|4.8e-6|1.1e-6|1.1e-6#, plastic_strain: 4.8e-6
    set Grain size                                = 1 #1e-3
    set Grain size exponents for diffusion creep  = 0   #default

###    # Plasticity parameters
    set Angles of internal friction = background: 0|0|0|0#, plastic_strain:0            # Angle Mohr
    set Cohesions                   = background: 40e6|1e30|1e30|1e30#, plastic_strain:40e6

    #set Strain weakening mechanism = plastic weakening with plastic strain only
    #set Start plasticity strain weakening intervals = 1.e30
    #set Cohesion strain weakening factors = 1
    #set Friction strain weakening factors = 1
    #set End plasticity strain weakening intervals = 2.e30

  end
################################
  #subsection Composition reaction model

    # Above this depth the compositional fields react: The first field gets
    # converted to the second field. Units: $m$.
   # set Reaction depth                                 = 660.e3

    # Reference density $\rho_0$. Units: $kg/m^3$.
    #set Reference density                              = 3300

    # The value of the specific heat $C_p$. Units: $J/kg/K$.
    #set Reference specific heat                        = 1250

    # The reference temperature $T_0$. Units: $K$.
    #set Reference temperature                          = 293

    # The value of the thermal conductivity $k$. Units: $W/m/K$.
    #set Thermal conductivity                           = 4.7

    # The value of the thermal expansion coefficient $\beta$. Units: $1/K$.
    #set Thermal expansion coefficient                  = 3.0e-5 #Default: 2e-5

    # The temperature dependence of viscosity. Dimensionless exponent.
    #set Thermal viscosity exponent                     = 0.5

    # The value of the constant viscosity. Units: $kg/m/s$.
    #set Viscosity                                      = 1e22 #1e22 #Default: 5e24
    #The value of the width
    #set Reaction width                                 = 5e3 #double

    #The value of the smoothing
 #   set Smooth reaction                                = true #bool
 # end  

end

#######################
# Gravity model
subsection Gravity model
  set Model name = ascii data
  #set Model name = radial constant
  #subsection Radial constant
  #  set Magnitude = 10
  #end
#  set Model name = ascii data
#  subsection Ascii data model
#    set Data directory = /home/fregehrk/opt/aspect/data/gravity-model/prem.txt
#  end
end

subsection Formulation
  set Formulation = anelastic liquid approximation
end

#######################
# Postprocessing
subsection Postprocess
  set List of postprocessors = velocity statistics, basic statistics, temperature statistics, \
                               visualization, topography, composition statistics, depth average, heat flux statistics, velocity boundary statistics
  subsection Visualization
    set List of output variables = strain rate, stress, nonadiabatic temperature, depth, shear stress, vertical heat flux, heat flux map, material properties, named additional outputs
    set Output format = vtu
    set Time between graphical output = 1e7
    set Interpolate output = true
    set Number of grouped files = 0
    subsection Heat flux map
      set Output point wise heat flux = false  # default=false
    end
  end
  subsection Depth average
    set Number of zones = 10
    set Output format = txt
    set Time between graphical output = 1e7
    set Depth boundaries of zones = 0, 10000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 100000
  end
end

[tjhei]

I don't quite understand the problem. Can you describe exactly what you do?

Are you stopping and restarting the computation after a certain amount of time? It looks like you run from 900 to 966 or 967 and the code aborts before you restart again at 900. Can you increase the walltime? Or checkpoint every 50 timesteps?

[elodie-kendall]

Hi Timo. I am just running the model once- it is running to time step 966/967 and then I see in the log file it resets to time step 900. It stays as the same job name, and doesn't cancel/fail, it just keeps doing this loop

[tjhei]

There is no place in the code that jumps from step 967 to step 900. Your output also shows that the code is resuming from a snapshot. That happens when you run aspect again. To me it looks like you are running several times and the code crashes/hangs/ runs out of wall time on between.

[elodie-kendall]

Sorry my mistake, you're right Timo it was just incredibly slow and not checkpointing enough before a timeout. Thanks!