# A description of convection in a 3d spherical shell with 
# a prescribed initial condition based on the shear wave
# velocity model S20RTS.

# Define the number of space dimensions we would like to 
# work in:
set Dimension                              = 3 

# Specify the time you want to let the model run for in 
# years and the output directory. Here we only calculate
# the instantaneous solution.
set End time                               = 0
set Use years in output instead of seconds = true
set Output directory                       = output-S20RTS

# The following variables describe how the pressure should
# be normalized. Here, we choose a zero average pressure
# at the surface of the domain. The 'Surface pressure' and
# 'Adiabatic surface temperature' are used to compute
# the adiabatic reference profile.
set Pressure normalization                 = surface
set Surface pressure                       = 0
set Adiabatic surface temperature          = 1600

# Here we specify the residual tolerance for the linear solver.
subsection Solver parameters
  subsection Stokes solver parameters
    set Linear solver tolerance = 1e-4
  end
end

# Here we specify the geometry of the domain, which is 
# a spherical shell with inner radius of 3481km and 
# outer radius of 6371km
subsection Geometry model
  set Model name = spherical shell
  subsection Spherical shell
    set Inner radius  = 3481000
    set Outer radius  = 6371000
  end
end

# This section specifies the temperature at the boundary of 
# the domain. Here we set the temperature to be constant,
# but different from the reference temperature to approximate
# boundary layers.
subsection Boundary temperature model
  set Fixed temperature boundary indicators = top, bottom
  set List of model names = spherical constant
  subsection Spherical constant
    set Inner temperature = 2000
    set Outer temperature = 1000
  end
end

# This section describes the gravity field, which is pointing
# towards the Earth's center with the same magnitude of 10 m/s^2
# everywhere
subsection Gravity model
  set Model name = radial constant
  subsection Radial constant
    set Magnitude = 10
  end
end

# This section prescribes the initial condition in the temperature
# field, which is chosen as a scaled version of the S20RTS shear
# wave velocity model (Ritsema et al., 2000). S20RTS is defined
# by spherical harmonics up to degree 20 that are radially interpolated
# with a cubic spline. 
subsection Initial temperature model
  set Model name = S40RTS perturbation
    subsection S40RTS perturbation

# The two input options here are S20RTS or the higher resolution
# S40RTS (Ritsema et al., 2011). One can choose to remove the 
# degree 0 from these files so that the depth average value
# is zero. 
    set Initial condition file name       = S20RTS.sph
    set Remove degree 0 from perturbation = false

# The following parameters determine the scaling from shear wave 
# velocity perturbation to temperature differences. We chose the
# scaling to density perturbation as 0.15
    set Vs to density scaling             = 0.15
    set Thermal expansion coefficient in initial temperature scaling = 3e-5

# This specifies the background temperature to which we add the 
# temperature difference.
    set Reference temperature             = 1600
  end
end

# The material model is based on the simple material model, which assumes
# a constant density, and other parameters as stated below.
subsection Material model
  set Model name = simple
    subsection Simple model
    set Reference density                 = 3300
    set Viscosity                         = 1e21
    set Thermal expansion coefficient     = 3e-5
    set Reference temperature             = 1600
    set Thermal conductivity              = 4.125 
    set Reference specific heat           = 1250
  end
end

# For this calculation we only do 2 global refinement steps. This resolution
# is too low to fully resolve the mantle flow, however it does capture
# the main features.
subsection Mesh refinement
  set Initial global refinement          = 2
  set Initial adaptive refinement        = 0
end

# We assume free slip at the inner and outer boundary
subsection Boundary velocity model
  set Tangential velocity boundary indicators = top, bottom
end


# We output the density, velocity, dynamic topography, geoid and heat flux density
# for plotting. 
subsection Postprocess
  set List of postprocessors = geoid, velocity statistics, heat flux map, heat flux statistics, dynamic topography, visualization, basic statistics

  subsection Visualization
    set Output format                 = vtu
    set List of output variables      = geoid, dynamic topography, heat flux map, density,viscosity, gravity
    set Time between graphical output = 0
    set Number of grouped files       = 1 

# We only have dirichlet boundaries with tangential velocities, so we can
# increase the output resolution as described in the documentation of the 'heat
# flux map' postprocessor.
    subsection Heat flux map
      set Output point wise heat flux = true
    end
  end

  subsection Geoid
    set Also output the gravity anomaly  = true
  end
end