geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
223 stars 235 forks source link

How do I add a hyperdense layer around the CMB in the Steinberger model? #5147

Closed Francyrad closed 1 year ago

Francyrad commented 1 year ago

Dear users I edited the original steinberger.prm model to get a 2D spherical shell, and it works nice. Now i'd like to add an hyperdense layer around the CMB to simulate a possible origin of the Superplumes like shown in many studies that used Citcom.

This is the code that works, without an hyperdense layer:

set Dimension                              = 2
set End time                               = 3e8
set Adiabatic surface temperature          = 1613.0
set Output directory                       = output-steinberger

set Nonlinear solver scheme                = iterated Advection and Stokes
set Nonlinear solver tolerance             = 1e-3

# spherical shell

subsection Geometry model
  set Model name = spherical shell
  subsection Spherical shell
    set Inner radius  = 3480000
    set Outer radius  = 6370000
  end
end

  set Tangential velocity boundary indicators = top, bottom
end

subsection Nullspace removal
  set Remove nullspace = net rotation
end

# We use the gravity model from PREM.
subsection Gravity model
  set Model name = ascii data
end

# Set temperature
  set Fixed temperature boundary indicators = top, bottom
  set List of model names = spherical constant

  subsection Spherical constant
    set Inner temperature = 3773
    set Outer temperature = 273
  end
end

# The initial temperature model consists of an adiabatic profile,
# thermal boundary layers at the surface and the core-mantle boundary,
# and a small harmonic perturbation to initiate convection.
subsection Initial temperature model
  set List of model names = adiabatic, function

  subsection Function
    set Coordinate system   = spherical
    set Variable names      = r, phi
    set Function constants  = r0 = 3481000, r1 = 6371000
    set Function expression = 30 * sin((r-r0)/(r1-r0)*pi)*sin(16*phi) + 20 * sin(2*(r-r0)/(r1-r0)*pi)*sin(12*phi)
  end

  subsection Adiabatic
    set Age top boundary layer = 1e8
    set Age bottom boundary layer = 1e8
  end
end

subsection Material model
  set Model name = Steinberger

  subsection Steinberger model

    set Data directory                                = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
    set Lateral viscosity file name                   = temp-viscosity-prefactor.txt
    set Material file names                           = pyr-ringwood88.txt
    set Radial viscosity file name                    = radial-visc-simple.txt

    set Maximum lateral viscosity variation           = 1e3
    set Maximum viscosity                             = 5e23
    set Minimum viscosity                             = 1e20

    set Use lateral average temperature for viscosity = false

    set Latent heat                                   = false
  end
end

subsection Heating model
  set List of model names = adiabatic heating, shear heating
end

subsection Formulation
  set Mass conservation = projected density field
  set Temperature equation = real density
end

subsection Compositional fields
  set Number of fields                     = 1
  set Names of fields                     = density_field
  set Types of fields                     = density
  set Compositional field methods        = prescribed field
end

subsection Mesh refinement
  set Initial adaptive refinement        = 2
  set Initial global refinement          = 5
  set Refinement fraction                = 0.95
  set Coarsening fraction                = 0.05
  set Strategy                           = temperature, boundary, minimum refinement function
  set Time steps between mesh refinement = 500

  subsection Boundary
    set Boundary refinement indicators = top, bottom
  end

  subsection Minimum refinement function
    set Coordinate system   = depth
    set Variable names      = depth, phi
    set Function expression = if(depth<700000, 6, 5)
  end
end

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

  subsection Visualization
    set Output format                 = vtu
    set List of output variables      = material properties, nonadiabatic temperature, named additional outputs
    set Time between graphical output = 1e6
  end
end

subsection Solver parameters
  subsection Stokes solver parameters
    set GMRES solver restart length = 200
    set Number of cheap Stokes solver steps = 500
  end
end

subsection Checkpointing
  set Steps between checkpoint = 10
end

The difficulties that now i encounter is that i don’t know how to add a circular hyperdense compositional field layer around the CMB. In the user manual (pg 92) there is written that, to add an active compositional field, i’ve to add this to my script:

# This is the new part: We declare that there will
# be two compositional fields that will be
# advected along. Their initial conditions are given by
# a function that is one for the lowermost 0.2 height
# units of the domain and zero otherwise in the first case,
# and one in the top most 0.2 height units in the latter.
subsection Compositional fields
set Number of fields = 2
end
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = if(y<0.2, 1, 0) ; if(y>0.8, 1, 0)
#Numbers will be bigger if i use a box with real kilometers.
end
end

Anyway, in the Steinberger simulation that i posted there are already compositional fields lines inside:

subsection Compositional fields
  set Number of fields                     = 1
  set Names of fields                     = density_field
  set Types of fields                     = density
  set Compositional field methods        = prescribed field
end

So i don't really get how to add this circular compositional field around the CMB honestly... I really need help.

Thanks in advance Francesco

Francyrad commented 1 year ago

Solved using 2 perplex files, particles and prescribed fields to both