# This is a cookbook which simulates a simple 2D subduction system with
# a hydrated sediment layer, MORB layer, gabbro layer, and peridotite
# layer. The layers dehydrate as they are subducted, with the solid-fluid
# reactions governed by an approximation published by Tian et al., 2018.
set World builder file                         = 2D_tian_slab.json
set Adiabatic surface temperature              = 1600
set Nonlinear solver scheme                    = iterated Advection and Stokes
set Nonlinear solver failure strategy          = cut timestep size
# set Nonlinear solver failure strategy          = continue with next timestep
set Output directory                           = cut_timestep_size_test_no_cutting
set Max nonlinear iterations                   = 20
set Nonlinear solver tolerance                 = 1e-5
set Dimension                                  = 2
set End time                                   = 4000
set CFL number                                 = 0.1
set Surface pressure                           = 0
# If this is larger than the melt generation then that means that the fluid will instantly dehydrate
# Perhaps this is unstable, especially if the solver tolerance is 1? Make this 10x less than the
# reaction time step in the volatile material model.
set Maximum time step                          = 100
set Use operator splitting                     = true
set Resume computation                         = false

subsection Time stepping
  set List of model names = repeat on nonlinear solver failure
end

subsection Solver parameters
  subsection Stokes solver parameters
    set Stokes solver type                              = block GMG
    set GMRES solver restart length                     = 1000
    set Number of cheap Stokes solver steps             = 20000
    set Use full A block as preconditioner              = true
    set Linear solver tolerance                         = 1e-6
    set Linear solver A block tolerance                 = 1e-2
    set Maximum number of expensive Stokes solver steps = 0
  end
  subsection Operator splitting parameters
    set Reaction time step                     = 50
    set Reaction time steps per advection step = 50
  end
end

subsection Discretization
  subsection Stabilization parameters
    set beta  = 0.052 # 0.5
    set cR    = 0.11 # 1.0
  end
end

subsection Checkpointing
  set Steps between checkpoint = 100
end

subsection Compositional fields
  set Number of fields = 7
  set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone1
  set Compositional field methods = darcy field, field, field, field, field, field, static
end

subsection Temperature field
  set Temperature method = field
end

subsection Initial composition model
  set List of model names = world builder
  subsection World builder
    set List of relevant compositions = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone1
  end
end

subsection Boundary composition model
  set Fixed composition boundary indicators         = left, left lithosphere, top, right lithosphere, right, bottom
  set List of model names                           = initial composition
end

subsection Geometry model
  set Model name = box with lithosphere boundary indicators
  subsection Box with lithosphere boundary indicators
    set Lithospheric thickness = 200e3
    set X extent      = 1500e3
    set Y extent      = 600e3

    set X repetitions = 15
    set Y repetitions = 5
    set Y repetitions lithosphere = 2

    set Box origin X coordinate = 0
    set Box origin Y coordinate = 0
  end
end


subsection Gravity model
  set Model name = function
  subsection Function
    set Variable names      = x,y
    set Function expression = 0; -9.81
  end
end

subsection Formulation
  set Formulation = Boussinesq approximation 
end

subsection Initial temperature model
  set Model name = world builder
end

subsection Boundary temperature model
  set Fixed temperature boundary indicators = top, bottom
  set List of model names = initial temperature
end


subsection Material model

  set Model name = reactive fluid transport
  set Material averaging                                = harmonic average only viscosity

  subsection Reactive Fluid Transport Model
    set Base model                                       = visco plastic
    set Reference fluid density                          = 1000
    set Shear to bulk viscosity ratio                    = 0.1
    set Reference fluid viscosity                        = 1
    set Reference permeability                           = 1e-6
    set Exponential fluid weakening factor               = 30
    set Fluid compressibility                            = 0
    set Fluid reaction time scale for operator splitting = 5e4
    set Fluid-solid reaction scheme                      = tian approximation
    set Maximum weight percent water in peridotite       = 2
    set Maximum weight percent water in gabbro           = 1
    set Maximum weight percent water in MORB             = 2
    set Maximum weight percent water in sediment         = 3
  end

  subsection Visco Plastic
    set Viscosity averaging scheme                        = harmonic
    set Viscous flow law                                  = composite
    set Prefactors for diffusion creep                    = background:4.5e-15, \
                                                            porosity:4.5e-15, \
                                                            bound_fluid:4.5e-15, \
                                                            peridotite:4.5e-15, \
                                                            gabbro:4.5e-15, \
                                                            MORB:4.5e-15, \
                                                            sediment:4.5e-15, \
                                                            weak_zone1:4.5e-15

    set Stress exponents for diffusion creep              = 1.0

    set Activation volumes for diffusion creep            = background:8.2e-6, \
                                                            porosity:8.2e-6, \
                                                            bound_fluid:8.2e-6, \
                                                            peridotite:8.2e-6, \
                                                            gabbro:8.2e-6, \
                                                            MORB:8.2e-6, \
                                                            sediment:8.2e-6, \
                                                            weak_zone1:8.2e-6

    set Activation energies for diffusion creep           = background:375e3, \
                                                            porosity:375e3, \
                                                            bound_fluid:375e3, \
                                                            peridotite:375e3, \
                                                            gabbro:375e3, \
                                                            MORB:375e3, \
                                                            sediment:375e3, \
                                                            weak_zone1:375e3

    set Prefactors for dislocation creep                  = background:7.4e-15, \
                                                            porosity:7.4e-15, \
                                                            bound_fluid:7.4e-15, \
                                                            peridotite:7.4e-15, \
                                                            gabbro:7.4e-15, \
                                                            MORB:7.4e-15, \
                                                            sediment:7.4e-15, \
                                                            weak_zone1:1e-50

    set Stress exponents for dislocation creep            = 3.5

    set Activation volumes for dislocation creep          = background:14e-6, \
                                                            porosity:14e-6, \
                                                            bound_fluid:14e-6, \
                                                            peridotite:14e-6, \
                                                            gabbro:14e-6, \
                                                            MORB:14e-6, \
                                                            sediment:14e-6, \
                                                            weak_zone1:14e-6

    set Activation energies for dislocation creep         = background:530e3, \
                                                            porosity:530e3, \
                                                            bound_fluid:530e3, \
                                                            peridotite:530e3, \
                                                            gabbro:530e3, \
                                                            MORB:530e3, \
                                                            sediment:530e3, \
                                                            weak_zone1:530e3

    set Angles of internal friction                       = 0

    set Cohesions                                         = background:1e50, \
                                                            porosity:1e50, \
                                                            bound_fluid:1e50, \
                                                            peridotite:1e50, \
                                                            gabbro:1e50, \
                                                            MORB:1e50, \
                                                            sediment:1e50, \
                                                            weak_zone1:10

    set Minimum viscosity                                 = 5e19
    set Maximum viscosity                                 = 5e24
    set Thermal expansivities                             = 3.1e-5
    set Thermal diffusivities                             = 1e-6
  end
end

subsection Mesh refinement
  set Coarsening fraction                      = 0.05
  set Refinement fraction                      = 0.8
  set Initial adaptive refinement              = 4
  set Initial global refinement                = 2
  set Strategy                                 = isosurfaces, composition threshold, minimum refinement function
  set Time steps between mesh refinement       = 2
  subsection Isosurfaces
    set Isosurfaces = max, max, bound_fluid: 0.005|0.04
  end

  # minimum of 4 global refinements
  subsection Minimum refinement function
    set Function expression = 2
  end

  subsection Maximum refinement function
    set Function expression = 7
  end

  # refine where the porosity is bigger than 1e-6
  subsection Composition threshold
    set Compositional field thresholds = 1e-6, 10.0, 10.0, 10.0, 10.0, 10.0, 0.5
  end
end

subsection Boundary velocity model
  set Tangential velocity boundary indicators = top, right lithosphere, right
  set Prescribed velocity boundary indicators = left lithosphere: function
  subsection Function
    set Function expression = 0.05; 0.0
  end
end

subsection Boundary traction model
  set Prescribed traction boundary indicators = left:initial lithostatic pressure, bottom: initial lithostatic pressure
  subsection Initial lithostatic pressure
    set Representative point = 0, 0
  end
end

subsection Melt settings
  set Include melt transport                  = false
end

subsection Postprocess
  set List of postprocessors = visualization, composition statistics, velocity statistics
  subsection Visualization
    set List of output variables      = density, darcy velocity, viscosity, thermal expansivity, named additional outputs
    set Output format                 = vtu
    set Time between graphical output = 0
    set Interpolate output            = true
    set Number of grouped files       = 0
  end
end