latent-heat.prm

############### 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 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
  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      = material properties

    subsection Material properties
      set List of material properties = density
    end
  end
end