Help improving the model speed

Hello, I am new-ish to aspect and I am trying to set up a 2d box model that is fairly simple to start with, that I intend to add complexity to (for example, plastic yielding, layered viscosity) hence why I am using the viscoplastic material model and for now set most parameters to zero.

The model runs, and the results look fine with visual inspection, but it is quite slow even for a simple setting. For example, I ran it on 112 cores for 3 hours and it only ran to ~ 4e7 years (the initial plumes only reached about half way through the domain). I was wondering if people have suggestions on how to improve the speed?

####  Global parameters
set Dimension                              = 2
set Start time                             = 0
set End time                               = 5e8
set Resume computation                     = auto
set Use years instead of seconds           = true
set Nonlinear solver scheme                = single Advection, single Stokes
set Nonlinear solver tolerance             = 5e-2
set Max nonlinear iterations               = 50
set CFL number                             = 0.5
set Output directory                       = output/
set Timing output frequency                = 1
set Pressure normalization                 = surface
set Surface pressure                       = 0

#### Parameters describing the model
set Adiabatic surface temperature          = 1873.

subsection Solver parameters
  subsection Stokes solver parameters
    set Linear solver tolerance                 = 1.0e-4
    set Number of cheap Stokes solver steps     = 300
    set Stokes solver type                      = block AMG
  end
  subsection Newton solver parameters
    set Max Newton line search iterations = 5
    set Max pre-Newton nonlinear iterations = 10
    set Maximum linear Stokes solver tolerance = 0.9
    set Nonlinear Newton solver switch tolerance = 1e-5
    set SPD safety factor = 0.9
    set Stabilization preconditioner = SPD
    set Stabilization velocity block = SPD
    set Use Eisenstat Walker method for Picard iterations = false
    set Use Newton failsafe = false
    set Use Newton residual scaling method = false
  end
  set Composition solver tolerance = 1e-12
  set Temperature solver tolerance = 1e-12
end

subsection Geometry model
  set Model name = box
  subsection Box
    set X extent = 6000e3
    set Y extent = 2980e3
    set X periodic = true
    set X repetitions = 8
    set Y repetitions = 4
  end
end

subsection Initial temperature model
  set Model name = function
  subsection Function
    set Coordinate system = cartesian
    set Variable names = x,y
    # Conductive geotherm plus a small lateral perturbation.
    # y=0 at the bottom and y=2980 km at the top in the box geometry.
    set Function constants = pi=3.141592653589793
    set Function expression = 1836 +  cos((y-1427e3)/2890e3*pi)*sin((4*x/3000e3)*300)
  end
end

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

  subsection Box
    set Bottom temperature = 2573
    set Top temperature    = 273
    set Left temperature   = 273
    set Right temperature  = 273
  end
end

# Gravity model
subsection Gravity model
  set Model name = vertical

  subsection Vertical
    set Magnitude = 9.81
  end
end

subsection Compositional fields
   set Number of fields = 1
   set Names of fields = total_strain
end

subsection Initial composition model
  set Model name = function
  subsection Function
    set Coordinate system = cartesian
    set Variable names = x, y
    set Function expression = 1e-8
  end
end

subsection Boundary velocity model
  set Tangential velocity boundary indicators = bottom, top
end

subsection Material model
  set Model name = visco plastic

  subsection Visco Plastic
    # Three depth-defined layers in the background material:
    #   0-150 km   : lithosphere
    #   150-400 km : asthenosphere
    #   >400 km    : mantle interior
    # These are implemented as depth-defined phases with zero Clapeyron slope,
    # so they behave as fixed layers in Cartesian geometry.
    set Phase transition depths            = background:150e3|400e3, total_strain:150e3|400e3
    set Phase transition widths            = background:50e3|50e3, total_strain:50e3|50e3
    set Phase transition temperatures      = background:1873|1873, total_strain:1873|1873
    set Phase transition Clapeyron slopes  = background:0|0, total_strain:0|0

    set Reference temperature = 1873
    set Minimum strain rate   = 1.e-20
    set Minimum viscosity     = 1e17
    set Maximum viscosity     = 2.5e25

    # No explicit temperature/pressure dependence in the viscous branch.
    set Adiabat temperature gradient for viscosity = 0.0

    # Density-related parameters
    set Thermal diffusivities = 1e-6
    set Heat capacities       = background:750|750|750, total_strain:750|750|750
    set Densities             = background:3300|3300|3300, total_strain:3300|3300|3300
    set Thermal expansivities = background:3e-5|3e-5|3e-5, total_strain:3e-5|3e-5|3e-5

    # Keep visco-plastic rheology, but make the viscous branch effectively layered
    # and nearly constant by removing T and P dependence from diffusion creep.
    # The constant viscosity prefactors impose the target layer viscosities.
    set Viscous flow law = diffusion
    set Activation energies for diffusion creep  = background:0|0|0, total_strain:0|0|0
    set Activation volumes for diffusion creep   = background:0|0|0, total_strain:0|0|0
    set Grain size exponents for diffusion creep = background:0|0|0, total_strain:0|0|0
    set Stress exponents for diffusion creep     = 1
    set Prefactors for diffusion creep           = background:2.0e-22|5.0e-20|1.0e-21, total_strain:2.0e-22|5.0e-20|1.0e-21
    #set Constant viscosity prefactors = background:100|0.1|10, total_strain:100|0.1|10

    # Plasticity is retained, but made weakly active so the initial viscosity profile
    # is controlled by the layered viscous branch rather than by yielding.
    set Yield mechanism            = drucker
    set Cohesions                  = background:1e9|1e9|1e9, total_strain:1e9|1e9|1e9
    set Angles of internal friction = background:0|0|0, total_strain:0|0|0
    # set Angles of plastic dilation  = background:0|0|0, total_strain:0|0|0
    set Reference strain rate      = 1e-15
    set Maximum yield stress       = 1e12
    set Viscosity averaging scheme = harmonic

    set Use plastic damper = true
    set Plastic damper viscosity = 1e21

    set Start prefactor strain weakening intervals = 0.
    set End prefactor strain weakening intervals   = 1.0
    set Prefactor strain weakening factors         = 1.0
    set Start plasticity strain weakening intervals = 0.
    set End plasticity strain weakening intervals   = 1.0
    set Cohesion strain weakening factors           = 1.0
    set Friction strain weakening factors           = 1.0

    set Strain weakening mechanism = total strain
    set Strain healing mechanism = temperature dependent
    set Strain healing temperature dependent prefactor = 250.
    set Strain healing temperature dependent recovery rate = 1.e-7
  end
  set Material averaging = harmonic average
end

subsection Formulation
  set Formulation = Boussinesq approximation
end

subsection Mesh refinement
  set Initial adaptive refinement         = 4
  set Initial global refinement           = 5
  set Strategy                            = temperature #, viscosity
  set Refinement criteria scaling factors = 1 #,1
  set Refinement criteria merge operation = max
  set Time steps between mesh refinement  = 10
end

subsection Postprocess
  set List of postprocessors = visualization, velocity statistics, temperature statistics, pressure statistics, heat flux statistics, velocity boundary statistics, mobility statistics

  subsection Visualization
    set Time between graphical output = 1e5
    set Number of grouped files = 1
    set List of output variables = material properties, strain rate, stress, temperature anomaly, heat flux map, vertical heat flux, compositional vector

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

subsection Nullspace removal
  # In Cartesian geometry remove translation null modes, not rotational ones.
  set Remove nullspace = linear x momentum, linear y momentum
end

subsection Checkpointing
  set Time between checkpoint = 500
end


Thanks in advance!

Did you read through this section in the manual? Making ASPECT run faster — ASPECT 3.1.0-pre