Skip to content

Commit

Permalink
add new Setup
Browse files Browse the repository at this point in the history
  • Loading branch information
@tgeffroy2
tgeffroy2 committed Feb 19, 2024
1 parent 0e1a31c commit 868ee8b
Show file tree
Hide file tree
Showing 2 changed files with 346 additions and 0 deletions.
217 changes: 217 additions & 0 deletions SETS/RiverTomCrouteinf.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,217 @@
#include "math.h"
#include "mdoodz.h"
#include "stdbool.h"
#include "stdlib.h"

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

double SetSurfaceZCoord(MdoodzInput *instance, double x_coord) {
double TopoLevel = 0.0;
// if (fabs(x_coord) < instance->model.surf_Winc/2.0) TopoLevel = -500.0/instance->scaling.L;
return TopoLevel;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

int SetPhase(MdoodzInput *instance, Coordinates coordinates) {
const double H_crust = instance->model.user1 / instance->scaling.L;
const double H_uppercrust = instance->model.user2 / instance->scaling.L;
int phase = 2;

if (coordinates.z> -H_crust) phase = 1;
if (coordinates.z> -H_uppercrust) phase = 0;

return phase;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

int SetDualPhase(MdoodzInput *input, Coordinates coordinate, int phase) {

int il, dual_phase = phase;
double Lx = input->model.xmax - input->model.xmin;
double Lz = input->model.zmax - input->model.zmin;
double Ax, Az;

// Set checkerboard for phase 0
Ax = cos( 50.0*M_PI*coordinate.x / Lx );
Az = sin( 80.0*M_PI*coordinate.z / Lz );
if ( ( (Az<0.0 && Ax<0.0) || (Az>0.0 && Ax>0.0) ) && dual_phase==0 ) {
dual_phase += input->model.Nb_phases;
}

if ( ( (Az<0.0 && Ax<0.0) || (Az>0.0 && Ax>0.0) ) && dual_phase==1 ) {
dual_phase += input->model.Nb_phases;
}

return dual_phase;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

double SetTemperature(MdoodzInput *instance, Coordinates coordinates) {
const double lithosphereThickness = instance->model.user1 / instance->scaling.L;
const double surfaceTemperature = 273.15 / instance->scaling.T;
const double mantleTemperature = (instance->model.user0 + 273.15) / instance->scaling.T;
const double particleTemperature = ((mantleTemperature - surfaceTemperature) / lithosphereThickness) * (-coordinates.z) + surfaceTemperature;
if (particleTemperature > mantleTemperature) {
return mantleTemperature;
} else {
return particleTemperature;
}
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

double SetHorizontalVelocity(MdoodzInput *instance, Coordinates coordinates) {
return -coordinates.x * instance->model.bkg_strain_rate;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

double SetVerticalVelocity(MdoodzInput *instance, Coordinates coordinates) {
return coordinates.z * instance->model.bkg_strain_rate;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coordinates) {
SetBC bc;
const double Lx = instance->model.xmax - instance->model.xmin;
const double V_tot = Lx * instance->model.bkg_strain_rate;
const double x = coordinates.x, z = coordinates.z;

// Evaluate velocity of W and E boundaries
const double VxW = 0.5*V_tot;
const double VxE = -0.5*V_tot;

// Assign BC values
if (position == N || position == S || position == NW || position == SW || position == NE || position == SE) {
bc.type = constant_shear_stress;
bc.value = 0;
} else if (position == W) {
bc.type = constant_velocity;
bc.value = VxW;
} else if (position == E) {
bc.type = constant_velocity;
bc.value = VxE;
} else {
bc.type = inside;
bc.value = 0.0;
}
return bc;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coordinates) {
SetBC bc;
const double Lx = instance->model.xmax - instance->model.xmin;
double V_tot = Lx * instance->model.bkg_strain_rate; // |VxW| + |VxE|

// Compute bottom boundary velocity and account for erosion
const double surf_Winc = instance->model.surf_Winc;
const double surf_Vinc = instance->model.surf_Vinc;
const double Vz_corr = surf_Winc*surf_Vinc / Lx;
const double VzS = (-0.5 * V_tot * (instance->topo_height->west - instance->model.zmin) / Lx) -(0.5 * V_tot * (instance->topo_height->east - instance->model.zmin) / Lx) + Vz_corr;

// Set boundary nodes types and values
if (position == W || position == SW || position == NW ) {
bc.type = constant_shear_stress;
bc.value = 0.0;
} else if ( position == E || position == SE || position == NE) {
bc.type = constant_shear_stress;
bc.value = 0.0;
} else if (position == S || position == N) {
bc.type = constant_velocity;
bc.value = VzS;
} else {
bc.type = inside;
bc.value = 0.0;
}
return bc;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

SetBC SetBCT(MdoodzInput *instance, POSITION position, double particleTemperature) {
SetBC bc;
double surface_temperature = zeroC / instance->scaling.T;
double mantle_temperature = (1330. + zeroC) / instance->scaling.T;
if (position == S) {
bc.type = constant_temperature;
bc.value = mantle_temperature;
}
if (position == free_surface || position == N) {
bc.type = constant_temperature;
bc.value = surface_temperature;
}
if (position == W || position == E) {
bc.type = constant_heatflux;
bc.value = 0.;
}
return bc;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

void AddCrazyConductivity(MdoodzInput *input) {
int *asthenospherePhases = (int *) malloc(sizeof(int));
CrazyConductivity *crazyConductivity = (CrazyConductivity *) malloc(sizeof(CrazyConductivity));
asthenospherePhases[0] = 3;
crazyConductivity->phases = asthenospherePhases;
crazyConductivity->nPhases = 1;
crazyConductivity->multiplier = 1000;
input->crazyConductivity = crazyConductivity;
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

int main() {
MdoodzSetup setup = {
.BuildInitialTopography = &(BuildInitialTopography_ff){
.SetSurfaceZCoord = SetSurfaceZCoord,
},
.SetParticles = &(SetParticles_ff){
.SetPhase = SetPhase,
.SetTemperature = SetTemperature,
.SetHorizontalVelocity = SetHorizontalVelocity,
.SetVerticalVelocity = SetVerticalVelocity,
.SetDualPhase = SetDualPhase,
},
.SetBCs = &(SetBCs_ff){
.SetBCVx = SetBCVx, // manuel
.SetBCVz = SetBCVz,
.SetBCT = SetBCT,
},
.MutateInput = AddCrazyConductivity,

};
RunMDOODZ("RiverTomCrouteinf.txt", &setup);
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/
129 changes: 129 additions & 0 deletions SETS/RiverTomCrouteinf.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,129 @@
/**** RESTART ****/
istep = 00020
irestart = 0

/**** OUTPUT FILES ****/
writer = 1
writer_step = 50
writer_markers = 0

/**** SCALES ****/ Scaling parameters for
eta = 1.e+22 /viscosity
L = 1.e4 /length
V = 1.0e-9 /velocity
T = 1.e+02 /temperature

/**** SPACE****/
Nx = 251
Nz = 81
xmin = -250.000000e3
xmax = 250.000000e3
zmin = -150.000000e3
zmax = 10.000000e3

/**** TIME ****/
Nt = 100000
dt = 2.63e10
constant_dt = 0
Courant = 0.25
RK = 4

/**** SWITCHES ****/
mechanical = 1
pure_shear_ALE = -1 /1: box stretched at constant strain rate; -1 box of constant size
elastic = 1
thermal = 1
free_surface = 1
free_surface_stab = 1.0
initial_cooling = 1
subgrid_diffusion = 2
shear_heating = 1
adiab_heating = 0
track_T_P_x_z = 1 / initial P and T field record
conserv_interp = 0

/**** SURFACE PROCESSES ****/
surface_processes = 5 /0 = none; 1 = diffusion only; 2 = fills instantaneously the basin up to the base level; 3 = diffusion + source term
surf_diff = 6e-6 /topography diffusion coefficient
surf_Winc = 1e3
surf_Vinc = 5e-9

/**** SETUP DEPENDANT ****/
bkg_strain_rate = 7.5e-16 / Background strain rate [s-1]
user0 = 1330.0 / Mantle temperature
user1 = 60.0e3 / Crust thickness
user2 = 20.0e3 /Upper crust Thickness


/**** GRAVITY ****/
gx = 0.0000
gz = -9.81

/**** PHASE PROPERTIES ****/
Nb_phases = 3

ID = 0 / upper Crust
rho = 2800.00 / ref. density
G = 3e10 / shear modulus
alp = 30.0e-6 / thermal expansivity
bet = 1e-11 / compressibility
Cv = 1050 / heat capacity
k = 2.0 / thermal conductivity
Qr = 1e-7 / radiogenic heat production
C = 5e7 / cohesion
phi = 30.00 / friction angle
eta_vp = 1.333e21 / viscoplastic viscosity
pwlv = 10 / power law viscosity -> database flow_laws.c
density_model = 3

ID = 1 / Lower Crust
rho = 2800.00 / ref. density
G = 3e10 / shear modulus
alp = 30.0e-6 / thermal expansivity
bet = 1e-11 / compressibility
Cv = 1050 / heat capacity
k = 2.0 / thermal conductivity
Qr = 1e-7 / radiogenic heat production
C = 5e7 / cohesion
phi = 30.00 / friction angle
eta_vp = 1.333e21 / viscoplastic viscosity
pwlv = 10 / power law viscosity -> database flow_laws.c
density_model = 3

ID = 2 / Mantle
rho = 3260.00 / ref. density
G = 3e10 / shear modulus
alp = 30.0e-6 / thermal expansivity
bet = 1e-11 / compressibility
Cv = 1050 / heat capacity
k = 3.0 / thermal conductivity
Qr = 1.0e-10 / radiogenic heat production
C = 5e7 / cohesion
phi = 30.00 / friction angle
eta_vp = 1.333e21 / viscoplastic viscosity
pwlv = 40 / power law viscosity -> database flow_laws.c
linv = 40 / linear creep viscosity -> database flow_laws.c
expv = 40 / exponential creep viscosity -> database flow_laws.c
density_model = 3

/**** PARTICLES ****/
Nx_part = 4
Nz_part = 4
min_part_cell = 16

/**** LINEAR SOLVER ****/
lin_abs_mom = 5.0e-8
lin_abs_div = 5.0e-8
lin_rel_mom = 5.0e-8
lin_rel_div = 5.0e-8
penalty = 1e3

/**** NON-LINEAR ITERATIONS ****/
Newton = 1
Newton2Picard = 1
nit_max = 20 // for high resolution use from 10 to 20
line_search = 0
nonlin_abs_mom = 5.0e-8
nonlin_abs_div = 5.0e-8
min_eta = 5.0e18
max_eta = 1.0e25

0 comments on commit 868ee8b

Please sign in to comment.