############### Global parameters
# We use a 2d setup. Since it takes some time for
# the model to reach a steady state we set the 
# end time to approximately 15 billion years. 

set Dimension = 2

set Start time                             = 0
set End time                               = 5e17
set Use years in output instead of seconds = false

set Output directory                       = output-latent-heat


subsection Geometry model
  set Model name = box

  subsection Box
    set X extent = 1000000
    set Y extent = 1000000
  end
end


subsection Gravity model
  set Model name = vertical
  subsection Vertical
    set Magnitude = 10.0
  end
end


subsection Heating model
  # As we only want to look at the effects of latent heating, we disable all
  # the other heating terms. 
  set List of model names =  latent heat
end


############### Boundary conditions
# We only fix the temperature at the upper boundary, the other boundaries
# are isolating.
subsection Boundary temperature model
  set Fixed temperature boundary indicators = top
  set List of model names = box

  subsection Box
    set Top temperature = 1000
  end
end

# To guarantuee a steady downward flow, we fix the velocity
# at the top and bottom, and set it to free slip on the sides. 
subsection Boundary velocity model
  set Prescribed velocity boundary indicators = bottom:function, top:function
  set Tangential velocity boundary indicators = left, right

  subsection Function
    set Function expression = 0;-2.1422e-11
    set Variable names      = x,y
  end
end

subsection Initial temperature model
  set Model name = function
  subsection Function
    set Function expression = 1.0e3
    set Variable names      = x,y
  end
end


subsection Material model
  set Model name = latent heat
  subsection Latent heat

    # The change of density across the phase transition. Together with the
    # Clapeyron slope, this is what determines the entropy change.
    set Phase transition density jumps                 = 115.6
    set Corresponding phase for density jump           = 0

    # If the temperature is equal to the phase transition temperature, the 
    # phase transition will occur at the phase transition depth. However, 
    # if the temperature deviates from this value, the Clapeyron slope 
    # determines how much the pressure (and depth) of the phase boundary
    # changes. Here, the phase transition will be in the middle of the box
    # for T=T1. 
    set Phase transition depths                        = 500000
    set Phase transition temperatures                  = 1000
    set Phase transition Clapeyron slopes              = 1e7

    # We set the width of the phase transition to 5 km. You may want to 
    # change this parameter to see how latent heating depends on the width
    # of the phase transition. 
    set Phase transition widths                        = 5000

    set Reference density                              = 3400
    set Reference specific heat                        = 1000
    set Reference temperature                          = 1000
    set Thermal conductivity                           = 2.38

    # We set the thermal expansion amd the compressibility to zero, so that 
    # all temperature (and density) changes are caused by advection, diffusion 
    # and latent heating. 
    set Thermal expansion coefficient                  = 0.0
    set Compressibility                                = 0.0

    # Viscosity is constant. 
    set Thermal viscosity exponent                     = 0.0
    set Viscosity                                      = 8.44e21
    set Viscosity prefactors                           = 1.0, 1.0
    set Composition viscosity prefactor                = 1.0
  end
end


subsection Mesh refinement
  set Initial adaptive refinement        = 0 
  set Initial global refinement          = 7 
  set Time steps between mesh refinement = 0
end


subsection Discretization
  subsection Stabilization parameters
    # The exponent $\alpha$ in the entropy viscosity stabilization. Units:
    # None.
    set alpha = 2

    # The $\beta$ factor in the artificial viscosity stabilization. An
    # appropriate value for 2d is 0.052 and 0.078 for 3d. Units: None.
    set beta  = 0.078

    # The $c_R$ factor in the entropy viscosity stabilization. Units: None.
    set cR    = 0.5   # default: 0.11
  end
end


subsection Postprocess
  set List of postprocessors = visualization

  subsection Visualization
    set Number of grouped files       = 0
    set Output format                 = vtu

    # We are only interested in the last timestep (when the system hast reached
    # a steady state). For following the development of the system or checking
    # if the solution already reached steady state, this parameter can be set 
    # to a smaller value.
    set Time between graphical output = 5e17
    set List of output variables      = density
  end
end