!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 8888888b. .d88888b. 888 888 d8888 8888888b. |
! 888 "Y88b d88P" "Y88b 888 888 d88888 888 Y88b |
! 888 888 888 888 888 888 d88P888 888 888 |
! 888 888 888 888 888 888 d88P 888 888 d88P |
! 888 888 888 888 888 888 d88P 888 8888888P" |
! 888 888 888 888 888 888 d88P 888 888 T88b |
! 888 .d88P Y88b. .d88P Y88b. .d88P d8888888888 888 T88b |
! 8888888P" "Y88888P" "Y88888P" d88P 888 888 T88b |
! |
! Input file for DOUAR-WSMP version 0.2 - 23/03/2011 |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
This file is read by two subroutines in the code:
- read_controlling_parameters
- read_input_file
The indications between [] indicate the type of the read parameter. It can be
an integer [int], a real*8 [dp], a character chain [char] or a boolean [bool].
CONTROLLING PARAMETERS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]debug is a parameter that switches on/off various prints and outputs
(the level of printing for error, warning and messages of the solver id%cntl(4)
is set to the debug value).
if debug is zero, no debugging
if debug is equal to one, this triggers the terminal display of some key
parameters
if debug is equal to two, same as debug equal to one, and the code produces
various vtk files in the DEBUG subdirectories (surfaces, olsf, ...), as well as
cross sections for all nonlinear iterations. Careful, memory consuming!
debug = 1
[bool]doDoRuRe is a flag that triggers the production of output files needed to
produce the DoRuRes. DoRuRe stands for 'Douar Run Report'.
doDoRuRe = F
[bool]compute_qp_gram triggers the production of qpgrams for every grid.
compute_qpgram = F
[bool]compute_reaction_forces toggles on/off the reaction forces computations.
compute_reaction_forces = F
RESTART
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]irestart is a restart flag; if irestart is not 0, the run will restart from
an output file given by [char]restartfile and at step.
irestart = 0
restartfile = OUT/time_0000.bin
TIMESTEPPING
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[dp]dt is the time step length (if dt is negative, courant conidition is used
and automatic time stepping is turned on)
dt = 1.25d-3
[int]nstep is the number of time steps
nstep = 50
[dp]courant is only used when dt is negative; it determines the size of the time
step from the maximum value of the velocity field amplitude. The time step s
the product of courant by the ratio of the smallest leaf size by the maximum
velocity.
courant = .5d0
[bool]normaladvect is a flag used to determine which algorithm to use to
calculate the new geoletry of the normals to the surfaces at the nodes on the
surfaces
if normaladvect is T, the normals are advected using the velocity gradient
if normaladvect is F, the normals are re-computed from the geometry of the
surface
normaladvect = T
GRID ITERATIONS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]griditer is a flag that allows for nonlinear iterations; when positive, a
fixed number (griditer) of iterations is permitted; when negative, the number of
nonlinear iterations is determined by a convergence criterion.
griditer = -100
[dp]octree_refine_ratio is the threshold value used to determine whether the
octree has converged or not. the larger the value, the less stringent the test.
octree_refine_ratio = .025d0
NONLINEAR ITERATIONS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]nonlinear_iterations is the maximum number of nonlinear iterations (i.e.
the iterations on a given constant grid)
if nonlinear_iterations is positive, it simply is the number of nonlinear
iterations performed for each grid. When negative it indicates an upper bound
of nonlinear iterations, but the actual number of nonlinear iterations is
determined by a convergence criterion (see the 'tol' parameter)
nonlinear_iterations = -100
[int]nb_iter_nl_min is the number of nonlinear iterations required for velocity
convergence when using the convergence criterion tolerance value below. Note
that this is different than simply requiring a given number of nonlinear
iterations as the velocity field must also be converged within the tolerance
below.
nb_iter_nl_min = 0
[dp]tol is the relative tolerance used to estimate convergence on the computed
velocity field
tol = 0.01d0
[bool]adaptive_tol is a flag that toggles on/off the evolution of the tol
parameter with the grid level: when velocity convergence is reached on a grid,
the latest meaure of the velocity difference between the two last obtained
solutions is put in tol, so that on the following generated grid, the solution
reaches at least the same level of convergence. It allows to start with a not
too stringent value of tol at uniform octree level that evolves with every
grid, assuming that increasing the level of refinement of the octree allows to
better capture the solution, hence allowing a tighter convergence.
adaptive_tol = F
OCTREES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]leveluniform_oct is the level of uniform discretization of space; note that
a level is a power of two used to divide the unit cube
leveluniform_oct = 4
[int]levelmax_oct is maximum level of octree discretization
levelmax_oct = 5
[double]vex is the vertical exaggeration scaling factor that allows for
variable aspect ratio elements. The value should range between 0-1, which
will be multiplied by the vertical scale of the model domain. For example, a
value of 0.25 would correspond to a model of dimensions 1x1x0.25, where the
elements are 1/4 as tall as they are wide. NOTE: vex MUST be larger than the
highest surface/element in the model at all times.
vex = 1.d0
[bool]ismooth is a flag to impose an additional level of smoothing after
refinement for the surfaces and strain rate. It ensures that no leaf is flanked
by other leaves differing by more than 1 level of refinement
ismooth = F
[int]noctreemax is the maximum size of any octree used in all computations
noctreemax = 10000000
[dp]refine_ratio is used to determine octree refinement based on a given
criterion. All leaves where the criterion is larger than refine_ratio times the
maximum of this criterion are refined
refine_ratio = -5500
[int]refine_criterion determines which refinement algorithm is to be used.
Several criteria exist for the refinement of the osolve octree. 1 is the second
invariant of the deviatoric strain-rate tensor; 2 is the sum of the squares of
the diagonal terms of the deviatoric strain-rate tensor; 3 is the second
invariant of the deviatoric strain rate tensor timses the leaf size. any other
value sets the criterion to zero and leads to no refinement.
refine_criterion = 0
[int]initial_refine_level is the initial level at which the refinement of the
octree will be performed. it has to be smaller than levelmax_oct
this should be used (different from levelmax_oct) in case the flow is very
localized (nonlinear/plastic analysis)
initial_refine_level = 4
[bool]renumber_nodes is a flag that can toggle on/off the renumbering of the
nodes by mean of Sloan's algorithm (T/F)
renumber_nodes = T
BOUNDARY CONDITIONS AND VELOCITY CONSTRAINTS (much more coming in v0.2)
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool]2d_invariants forces the calculation of the strain rate tensor invariants
to two dimensions. This assumes no strain in the y direction, and sets eyy, exy
and eyz to zero.
invariants_2d = F
PRESSURE
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]smoothing_type is a parameter allows to choose which type of smoothing is
to be applied to the pressure field: 0 is none, 1 is center->nodes->center, 2 is
the same, but weighted by neighbouring elemental volumes, 3 is regular grid+SPH,
and 4 is SPH.
smoothing_type = 1
CLOUD
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]npmin and [int]npmax are used to update the 3D volumetric cloud. npmin
corresponds to the minimum number of particles in any leaf; npmax is the maximum
allowable number in any leaf
npmin = 1
npmax = 4
FEM + DIVFEM + MUMPS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]matrule determines the rule used to assign material properties to
elements that are cut by a surface. The default value (0) will use divFEM to
determine the element properties. A value of 1 will assign the properties of
the material that occupies the majority of the element, whereas a value of 2
will assign the properties of the volumetric minority of the element. NOTE:
divFEM is always used for any element containing void material.
matrule = 1
[int]levelcut is the number of levels used to estimate the volume integrals in
the divFEM algorithm for cut cells; by testing, it has been estimated that a
value of 2 is very accurate
levelcut = 2
[int]levelapprox is the number of levels used to estimate the remaining
integrals using an improved version of Marthijn's clever algorithm, usually 3 is
plenty
levelapprox = 3
[dp]penalty is a global penalty parameter used to impose the bad faces or
incompatible faces linear constraints
penalty = 1.d8
[bool]excl_vol is a parameter that toggles off the assumption that lsf's are
built on top of one another
excl_vol = F
TEMPERATURE
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[dp] ztemp is the height interval between which a linear temperature gradient is
set: the temperature is 1 at the bottom, and 0 at ztemp.
[dp] tempscale
calculate_temp = F
ztemp = 0.2d0
tempscale = 100.d0
ISOSTASY
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool]isostasy enables calculation of local isostasy with the equilibrium state
based on the initial geometry. Local isostasy will be used if this is true and
flexure is false.
isostasy = F
[bool]flexure enables calculation of flexural isostasy using the parameters
below. Note that either local or flexural isostasy will be used, based on this
value. Both must be true for flexure to be used.
flexure = F
[dp]elastic_plate_thickness is the thickness of the elastic plate used in the
flexural isostasy calculation, in meters (?)
elastic_plate_thickness = 20.d3
[dp]density_difference is the difference in density between the crust and mantle
in ??? I believe this is abs(rhoc-rhom)*g. (kg/m^3 * m/s^2)
density_difference = 3.d3
[bool]isobc enables modification of the velocity boundary conditions for certain
boundary condition types to account for changes in the Moho depth
isobc = F
MATERIALS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]nmat is number of materials
nmat = 1
[int]material0 determines what is the material above the 1st surface (the free
surface) if material0 is 0 then it is the void and the properties "0" are used
for that part of the model; otherwise the material is one of the material,
comprised between 1 and nmat
material0 = 0
[bool]bulkvisc determines whether the user will input a material penalty value
that scales with the shear viscosity or whether an independent bulk viscosity
will be listed. If true, the values listed for [dp]penaltyi below should be
bulk viscosity values for each material. Otherwise, the listed penalty value
will be multiplied by the shear viscosity to enforce the incompressibility
condition.
bulkvisc = F
[bool]init_e2d determines whether an initial e2d value will be used to put
plastic materials on yield in the first time step, or whether the standard
initial material viscosity will be used. If true, an initial e2d value will be
used to ensure plastic materials behave as desired.
init_e2d = T
[dp]e2d0 is the initial strain rate value to use if init_e2d is true. This
value will be used in only the first nonlinear iteration of the first time
step and should be in scaled DOUAR units.
e2d0 = 1.d-2
[dp]densityi, [dp]viscosityi and [dp]penaltyi are the density, viscosity and
incompressibility used for material i; there should be nmat sets of material
properties; there should also be a nil material if material0 has been set to 0
[dp]expon is the nonlinear viscosity exponent
[dp]diffusivity is the heat diffusivity
[dp]expansion is the coefficient of thermal expansion
[dp]heat is the heat production
[char]plasticity_type is the type of plasticity
- No -> no plasticity, purely (nonlinear) viscous material
- vM -> von Mises yield criterion
-> [dp]plasticity_1st_param is the yield strength
- DPI-VII-> Drucker-Prager type of yield criterion
DPI -> the yield locus passes through the outer apices of the Mohr-Coulomb
hexagonal cone
DPII -> the yield locus passes through the inner apices of the Mohr-Coulomb
hexagonal cone
DPIII -> The average case, intermediate to DPI and DPII (USE ONLY FOR TESTING;
this case can not easily be derived as a simplification of the
Mohr-Coulomb criterion)
DPIV -> Plane strain Drucker-Prager formulation for associative materials
(USE ONLY FOR TESTING)
-> [dp]plasticity_1st_param is the angle phi
-> [dp]plasticity_2nd_param is the cohesion c
-> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
which phi is modified
-> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
which phi is set to the value in plasticity_5th_param
-> [dp]plasticity_5th_param is the final phi value, if modified
- DPV -> Drucker-Prager type of yield criterion that does not calculate the
values of alpha and k from phi and c. They are specified here. Note
that strain softening cannot be used for this DP type.
-> [dp]plasticity_1st_param is alpha
-> [dp]plasticity_2nd_param is k
- DPVI -> Drucker-Prager type of yield criterion that uses the equation below
sigma_y = pressure * tan(phi) + c
(USE ONLY FOR TESTING)
-> [dp]plasticity_1st_param is the angle phi
-> [dp]plasticity_2nd_param is the cohesion c
-> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
which phi is modified
-> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
which phi is set to the value in plasticity_5th_param
-> [dp]plasticity_5th_param is the final phi value, if modified
- DPVII -> Drucker-Prager type of yield criterion that uses the same 2D
formulation that is used in Sopale
sigma_y = pressure*sin(phi) + c*cos(phi)
-> [dp]plasticity_1st_param is the angle phi
-> [dp]plasticity_2nd_param is the cohesion c
-> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
which phi is modified
-> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
which phi is set to the value in plasticity_5th_param
-> [dp]plasticity_5th_param is the final phi value, if modified
- MC -> Mohr-Coulomb type of yield criterion
-> [dp]plasticity_1st_param is the angle phi
-> [dp]plasticity_2nd_param is the cohesion c
-> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
which phi is modified
-> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
which phi is set to the value in plasticity_5th_param
-> [dp]plasticity_5th_param is the final phi value, if modified
[dp]ztrans is the elevation below which the material transition from the given
material to material transnum
[int]transnum is the material number used when z is below ztrans
density0 = 0.d0
viscosity0 = 1.d-5
penalty0 = 1.d8
expon0 = 1.d0
diffusivity0 = 1.d0
heat0 = 0.d0
activationenergy0 = 0.d0
plasticity_type0 = No
plasticity_1st_param0 = 0.d0
plasticity_2nd_param0 = 0.d0
plasticity_3rd_param0 = 0.5d0
plasticity_4th_param0 = 1.5d0
plasticity_5th_param0 = 0.d0
ztrans0 = -1.d0
transnum0 = 0
density1 = -1.d0
viscosity1 = 5.d4
penalty1 = 1.d8
expon1 = 1.d0
diffusivity1 = 2.d-3
heat1 = 0.d0
activationenergy1 = 0.d0
plasticity_type1 = MC
plasticity_1st_param1 = 15.d0
plasticity_2nd_param1 = 5.d-4
plasticity_3rd_param1 = 0.5d0
plasticity_4th_param1 = 1.5d0
plasticity_5th_param1 = 15.d0
ztrans1 = -1.d0
transnum1 = 1
[dp]viscositymin and viscositymax are bounds on the viscosity (if negative bound
is not imposed) These bounds are introduced to prevent the viscosity to reach
unrealistic values, especially when using non-linear (power-law or brittle)
rheologies
viscositymin = 5.d-8
viscositymax = 5.d9
SURFACES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]ns is number of surfaces to track
ns = 1
for each surface, one needs to define a levelt, itype, material, surface_type,
activation_time and surface_params.
- [int]levelt is the inital level for the particles on the surface; to be
accurate and avoid holes in the surface during definition of the lsf, one
should use levelt eq to levelmax_oct+1 for all surfaces as a minimum value;
- [int]itype should be 1 for foldable surfaces or 0 for nonfoldable surfaces;
- [int]material is the material type refering to the table of material available
(max nmat);
- [dp]activation_time is the time the surface becomes active (before that time,
it is glued to the 0th surface). this parameter is useful when defining
stratigraphic horizons; default is -1, ie the surface is not glued to the free
surface
- [int]surface_type is comprised between 1 and 13.
1 corresponds to a flat surface,
-> surface_param_01 is the z level
2 to a rectangular emboss,
-> surface_param_01 is the z level
-> surface_param_02 and 03 are x1,x2
-> surface_param_04 and 05 are y1,y2
-> surface_param_06 is the thickness
3 to a convex spherical emboss,
-> surface_param_01 is the z level
-> surface_param_02 and 03 are x0,y0
-> surface_param_04 is the radius
4 to concave spherical emboss,
-> surface_param_01 is the z level
-> surface_param_02 and 03 are x0,y0
-> surface_param_04 is the radius
5 to a double rectangular emboss,
-> surface_param_01 is the z level
-> surface_param_02 and 03 are x1,x2
-> surface_param_04 and 05 are x3,x4
-> surface_param_06 and 07 are y1,y2
-> surface_param_08 and 09 are y3,y4
-> surface_param_10 is the thickness
6 to a sinus,
-> surface_param_01 is the z level
-> surface_param_02 is the wavelength
-> surface_param_03 is the amplitude
7 to a noisy surface,
-> surface_param_01 is the z level
-> surface_param_02 is the noise amplitude
8 to a double sinus,
-> surface_param_01 is the z level
-> surface_param_02 is the x-wavelength
-> surface_param_03 is the x-amplitude
-> surface_param_04 is the y-wavelength
-> surface_param_05 is the y-amplitude
9 to a cosinus,
-> surface_param_01 is the z level
-> surface_param_02 is the wavelength
-> surface_param_03 is the amplitude
10 to a 2D embankment,
-> surface_param_01 is z0
-> surface_param_02 is y0
-> surface_param_03 is psi
-> surface_param_04 is the thickness
11 to a 2D hill,
-> surface_param_01 is z level
-> surface_param_02 is the height
-> surface_param_03 is y0
-> surface_param_04 is width
12 to a rectangular emboss with specified slope,
-> surface_param_01 is the z level
-> surface_param_02 is x1
-> surface_param_03 is x2
-> surface_param_04 is y1
-> surface_param_05 is y2
-> surface_param_06 is the thickness
-> surface_param_07 is the slope
13 to a 2D embankment of finite length and with 2 kinks
-> surface_param_01 is the z level
-> surface_param_02 is the thickness (zthick)
-> surface_param_03 is slope angle (psi)
-> surface_param_04 is x value for the center of the slope (x0)
-> surface_param_05 is the x width of the narrow end of the plateau (y0w)
-> surface_param_06 is the y position for the start of the slope (y0)
-> surface_param_07 is the y position of the first kink (y1)
-> surface_param_08 is the angle of the first kink (theta1)
-> surface_param_09 is the y position of the second kink (y2)
-> surface_param_10 is the angle of the second kink (theta2)
- [int]leveloct is the level at which the octree will be refined in the vicinity
of the surface.
levelt1 = 6
itype1 = 1
surface_type_1 = 1
rand1 = F
surface_param_01_1 = 0.2d0
material1 = 1
activation_time_1 = -1.
leveloct1 = 5
stretch1 = 1.5d0
anglemax1 = 180.d0
criterion1 = 1
anglemaxoctree1 = 180.d0
spread_surface_points1 = 1
[int]niter_move is the number of iterations used to update particle positions
using an implicit, mid-point algorithm (default is 10)
niter_move = 10
[dp]stretch is the maximum allowed increase in linear length between two
initially adjacent particles on any surface; when this stretch is achieved, a
new particle is inserted on the surface, half-way along the stretched edge
stretch = 1.5d0
[dp]anglemax is the maximum allowed angle between two adjacent normals
when the angle is reached a new point is inserted bewteen the two points to
reduce the angle between the two normals
anglemax = 180.d0
[int]criterion is criterion used to define the octree in the vicinity of the
sufaces; criterion 1 corresponds to imposing that all leaves that are cut by any
of the surfaces must be at level levelmax_oct; criterion 2 corresponds to
imposing that discretization is proportional to the curvature of the surface;
curvature is calculated from the local divergence of the normals. criterion 3
corresponds to imposing that all leaves that contain at least one particle of
any surface is at levelmax_oct;
criterion = 2
[dp]anglemaxoctree is only defined for criterion 2; t is the maximum allowable
angle between two adjacent normals; if the angle is greater than anglemaxoctree,
the local octree leaves are forced to be at level levelmax_oct; otherwise they
are proportionally larger (smaller levels) (default is 10)
anglemaxoctree = 10.d0
REFINEMENT IN BOXES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]nboxes is the number of boxes in which the user imposes a set level of
discretization
nboxes = 1
for each box we need to specify the two end corners and the level
the syntax is [dp]boxnx0, [dp]boxnx1, [dp]boxny0, [dp]boxny1, [dp]boxnz0,
[dp]boxnz1 and boxnlevel where n is the box number
box1x0 = 0.d0
box1x1 = 1.d0
box1y0 = 0.d0
box1y1 = 1.d0
box1z0 = 0.d0
box1z1 = 0.25d0
box1level = 5
REFINEMENT ON CUBE FACES
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool]ref_on_faces toggles on/off the user imposed refinement on faces.
For each of the six faces of the cube, on defines the level at which the
desired area is to be refined. this area is given by bottom, top, left and
right coordinates that are comprised between 0 and 1
faces 1,2,3,4,5,6 respectively correspond to x=0,x=1,y=0,y=1,z=0,z=1
ref_on_faces = F
level_face1 = 5
b1 = .02
t1 = .51
l1 = .11
r1 = .81
level_face2 = 5
b2 = .45
t2 = .55
l2 = .46
r2 = .56
level_face3 = 5
b3 = .0
t3 = .2
l3 = .3
r3 = .7
level_face4 = 5
b4 = .4
t4 = .5
l4 = .14
r4 = .4
level_face5 = 6
b5 = 0.
t5 = 1.
l5 = 0.
r5 = 0.75
level_face6 = 5
b6 = .26
t6 = .56
l6 = .16
r6 = .86
EROSION
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool]erosion is a flag that toggles on/off the erosion.
erosion = T
if erosion is on, one also needs to set the erosion level/height, [dp]zerosion
Note that this is a first attempt at erosion; in future versions, DOUAR
should be easily linked to a surface processes model like CASCADE
zerosion = 0.2d0
If erosion is on, one also needs to define a length scale and a velocity scale
to properly translate the uplift rate produced by DOUAR into something that is
adequate for CASCADE;
[dp]length_scale is the scale of the unit model in km
[dp]velocity_scale is the velocity scale in km/Myr.
if length_scale is negative, erosion is assumed to be perfect (no call to
CASCADE is needed)
length_scale = 640.d0
velocity_scale = 50.d0
One then needs to define the erosion constants in CASCADE
[dp]fluvial_erosion is the fluvial erosion constant in 1/m^2
[dp]diffusion_erosion is the diffusion erosion constant in m^2/yr
(4d-2 32d-2)
fluvial_erosion = 4.d-2
diffusion_erosion = 32.d-2
One finally needs to specify the boundary conditions for CASCADE
if [int]baselevelx0 is set to 1 then the boundary at x=0 is set at baselevel
(water and sediment exit)
if [int]baselevelx1 is set to 1 then the boundary at x=1 is set at baselevel
(water and sediment exit)
if [int]baselevely0 is set to 1 then the boundary at y=0 is set at baselevel
(water and sediment exit)
if [int]baselevely1 is set to 1 then the boundary at y=1 is set at baselevel
(water and sediment exit)
baselevelx0 = 1
baselevelx1 = 1
baselevely0 = 0
baselevely1 = 1
MATRIX VISUALISATION
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[bool]visualise_matrix allows the user to turn on the visual representation of
the matrices used in the code. Be careful, the generated postscript files are
huge!
visualise_matrix = F
CROSS SECTIONS
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
[int]nsections is the number of cross-sections to be output. [dp]scale is the
scale used to produce the postscript file.
[int]xyz takes values 1, 2 or 3, and corresponds to planes defined by
x=constant, y=constant, and z=constant resp.
The [bool]flags are self explanatory. [char]colormap is the chosen colormap (jet
or hot). [int]ncolours is the number of colours used to produce the plot.
nsections = 0
xyz_1 = 1
slice_1 = 0.5001
flag_press_1 = T
flag_spress_1 = T
flag_e2d_1 = T
flag_e3d_1 = F
flag_strain_1 = F
flag_lode_1 = F
flag_crit_1 = F
flag_grid_1 = T
flag_mu_1 = F
flag_u_1 = F
flag_v_1 = F
flag_w_1 = F
flag_q_1 = F
flag_uvw_1 = F
flag_lsf_1 = F
flag_vfield_1 = F
flag_colour_1 = T
flag_plastic_1 = F
flag_temp_1 = F
flag_velvect_1 = T
scale_1 = 500.
colormap_1 = jet
ncolours_1 = 256
xyz_2 = 2
slice_2 = 0.901
flag_press_2 = T
flag_e2d_2 = T
flag_e3d_2 = F
flag_strain_1 = F
flag_lode_2 = F
flag_crit_2 = F
flag_grid_2 = F
flag_mu_2 = F
flag_u_2 = F
flag_v_2 = F
flag_w_2 = F
flag_q_2 = F
flag_uvw_2 = F
flag_lsf_2 = F
flag_vfield_2 = F
flag_colour_2 = T
flag_plastic_2 = F
flag_temp_1 = F
flag_velvect_2 = F
scale_2 = 800.
colormap_2 = jet
ncolours_2 = 256
xyz_3 = 3
slice_3 = 0.0010
flag_press_3 = F
flag_e2d_3 = F
flag_e3d_3 = F
flag_strain_1 = F
flag_lode_3 = F
flag_crit_3 = F
flag_grid_3 = F
flag_mu_3 = F
flag_u_3 = F
flag_v_3 = F
flag_w_3 = F
flag_q_3 = F
flag_uvw_3 = F
flag_lsf_3 = F
flag_vfield_3 = F
flag_colour_3 = T
flag_plastic_3 = F
flag_temp_1 = F
flag_velvect_3 = F
scale_3 = 400.
colormap_3 = jet
ncolours_3 = 256