Skip to content

Commit

Permalink
add new surface processes implementation (#125)
Browse files Browse the repository at this point in the history 
The current implementation uses a dm swarm to track the surface.
In this version, the only surface process available is a simple diffusion process.
Additionally, some parameters are no longer supported: `sp_dt`, `precipitation_profile_from_ascii`, and `climate_change_from_ascii`.
The output file for the surface level has been renamed to `surface_*.txt`.
  • Loading branch information
@rafaelmds
rafaelmds committed Sep 2, 2024
1 parent 6089f97 commit 055b4b3
Show file tree
Hide file tree
Showing 15 changed files with 1,687 additions and 1,627 deletions.
43 changes: 10 additions & 33 deletions src/header.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ PetscInt particles_per_ele = 81;
PetscReal theta_FSSA = 0.5;
PetscReal sub_division_time_step = 1.0;
PetscReal particles_perturb_factor = 0.5;
PetscInt sp_mode = 1;
PetscInt sp_mode = 0;
PetscReal Xi_min = 1.0E-14;
PetscReal random_initial_strain = 0;
PetscReal pressure_const = -1.0;
Expand All @@ -22,7 +22,6 @@ PetscScalar K_fluvial = 2.0E-7;
PetscScalar m_fluvial = 1.0;
PetscScalar sea_level = 0.0;
PetscScalar basal_heat = -1.0;
PetscReal sp_dt = 0.0;
PetscScalar sp_d_c = 0.0;
// Parameter file boolean variables
PetscInt WITH_NON_LINEAR = 0; // 1=True, 0=False
Expand All @@ -43,8 +42,6 @@ PetscInt bcv_extern = 0; // 1=True, 0=False
PetscInt binary_output = 0; // 1=True, 0=False
PetscInt sticky_blanket_air = 0; // 1=True, 0=False
PetscInt multi_velocity = 0; // 1=True, 0=False
PetscInt precipitation_profile = 0; // 1=True, 0=False
PetscInt climate_change = 0; // 1=True, 0=False
PetscInt free_surface_stab = 1; // 1=True, 0=False
PetscInt print_step_files = 1; // 1=True, 0=False
PetscInt RK4 = 0; // 0=Euler, 1=Runge-Kutta (not working yet!)
Expand All @@ -55,7 +52,6 @@ PetscInt high_kappa_in_asthenosphere = 0; // 1=True, 0=False
// Will be added to param.txt
PetscBool sp_surface_tracking = PETSC_FALSE; // PETSC_TRUE/PETSC_FALSE
PetscBool sp_surface_processes = PETSC_FALSE; // PETSC_TRUE/PETSC_FALSE
PetscBool set_sp_dt = PETSC_FALSE; // PETSC_TRUE/PETSC_FALSE
PetscBool set_sp_d_c = PETSC_FALSE; // PETSC_TRUE/PETSC_FALSE
PetscBool plot_sediment = PETSC_FALSE; // PETSC_TRUE/PETSC_FALSE
PetscBool a2l = PETSC_TRUE; // PETSC_TRUE/PETSC_FALSE
Expand Down Expand Up @@ -359,34 +355,6 @@ PetscInt cont_mv=0;

PetscInt print_step_aux;

Vec sp_surface_global;
Vec sp_surface_global_n;
Vec sp_surface_coords_global;
Vec sp_top_surface_global;
Vec sp_bot_surface_global;
Vec sp_mean_surface_global;
Vec sp_top_surface_global_n;
Vec sp_bot_surface_global_n;
PetscReal sp_eval_time;
PetscReal sp_last_eval_time;

long sp_n_profiles;
PetscScalar *topo_var_time;
PetscScalar *topo_var_rate;

Vec sp_surface_global_aux;

PetscScalar *global_surface_array_helper;
PetscScalar *global_surface_array_helper_aux;

PetscScalar prec_factor=1.0;

PetscScalar *var_climate_time;
PetscScalar *var_climate_scale;
PetscInt n_var_climate;
PetscInt cont_var_climate=0;


//Rescaled variables for non-dimensional scenarios
//(necessary to calculate the effective viscosity using
//dimensional pressure, temperature, strain rate and cummulative strain)
Expand All @@ -412,3 +380,12 @@ PetscReal pressure0_scaled;
PetscReal air_threshold_density;

PetscBool export_kappa = PETSC_FALSE;

// surface processes
DM dmcell_s;
DM dms_s;
PetscInt dms_s_ppe = 2; // dm swarm surface number of particles per element
PetscInt buffer_s;

// sp_mode
// 1 - diffusion
118 changes: 48 additions & 70 deletions src/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,15 +50,15 @@ PetscErrorCode calc_drho();
PetscErrorCode veloc_total(int dimensions);
PetscErrorCode rescaleVeloc(Vec Veloc_fut, double tempo);
PetscErrorCode multi_veloc_change(Vec Veloc_fut,double tempo);
PetscErrorCode sp_create_surface_vec();
PetscErrorCode sp_interpolate_surface_particles_to_vec();
PetscErrorCode evaluate_surface_processes();
PetscErrorCode sp_write_surface_vec(PetscInt i);
PetscErrorCode sp_destroy();
PetscErrorCode rescalePrecipitation(double tempo);
PetscErrorCode parse_options(int rank);
PetscErrorCode load_topo_var(int rank);

PetscErrorCode sp_create_surface_swarm_2d();
PetscErrorCode sp_move_surface_swarm(PetscInt dimensions, PetscReal dt);
PetscErrorCode sp_surface_swarm_interpolation();
PetscErrorCode sp_evaluate_surface_processes(PetscInt dimensions, PetscReal dt);
PetscErrorCode sp_update_surface_swarm_particles_properties();
PetscErrorCode sp_destroy();
PetscErrorCode sp_view_2d(DM dm, const char prefix[]);

int main(int argc,char **args)
{
Expand Down Expand Up @@ -96,34 +96,35 @@ int main(int argc,char **args)
seed = rank;

// Parse command line options
parse_options(rank);
ierr = parse_options(rank); CHKERRQ(ierr);

// Read ASCII files
reader(rank, "param.txt");
ierr = reader(rank, "param.txt"); CHKERRQ(ierr);

// Check if the number of interfaces in "param.txt" is higher than then number read from command line
if (seed_layer_set && seed_layer_size > (n_interfaces + 1)) {
PetscPrintf(PETSC_COMM_WORLD, "Error: The number of layers specified in command line \"-seed\" command is higher than the number in \"param.txt\".\n");
}

// Check surface processes configuration
if (sp_surface_processes == PETSC_TRUE && (sp_mode != 1)) {
PetscPrintf(PETSC_COMM_WORLD, "Error: Invalid \"sp_mode\" in \"param.txt\".\n");
exit(1);
}

PetscPrintf(PETSC_COMM_WORLD, "Number of seed layers: %d\n", seed_layer_size);
for (int k = 0; k < seed_layer_size; k++) {
PetscPrintf(PETSC_COMM_WORLD, "seed layer: %d - strain: %lf\n", seed_layer[k], strain_seed_layer[k]);
}
PetscPrintf(PETSC_COMM_WORLD, "\n");

if (sp_surface_processes && sp_surface_tracking && sp_mode == 1) {
load_topo_var(rank);
}

if (sp_mode == 2 && PETSC_FALSE == set_sp_d_c) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "-sp_mode 2 (diffusion) using default value: sp_d_c %e\n", sp_d_c); CHKERRQ(ierr);
} else if (sp_mode == 2) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "-sp_mode 2 (diffusion) using custom value: sp_d_c %e\n", sp_d_c); CHKERRQ(ierr);
} else if (sp_mode == 3) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "-sp_mode 3 (fluvial erosion) using K_fluvial: %e and sea_level %e\n", K_fluvial, sea_level); CHKERRQ(ierr);
} else if (sp_mode == 4) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "-sp_mode 4 (fluvial erosion mode 2) using K_fluvial: %e and sea_level %e\n", K_fluvial, sea_level); CHKERRQ(ierr);
// surface processes
if (sp_surface_processes == PETSC_TRUE && sp_mode == 1) {
if (set_sp_d_c == PETSC_FALSE) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "'sp_mode=1' (diffusion) using default value: sp_d_c %e\n", sp_d_c); CHKERRQ(ierr);
} else {
ierr = PetscPrintf(PETSC_COMM_WORLD, "'sp_mode=1' (diffusion) using custom value: sp_d_c %e\n", sp_d_c); CHKERRQ(ierr);
}
}

// Update elements aux constants
Expand Down Expand Up @@ -164,17 +165,17 @@ int main(int argc,char **args)
PetscPrintf(PETSC_COMM_WORLD,"Swarm: done\n");
}

// PetscPrintf(PETSC_COMM_SELF,"********** <rank:%d> <particles_per_ele:%d>\n", rank, particles_per_ele); -> conversar com Victor sobre
// PetscPrintf(PETSC_COMM_SELF,"********** <rank:%d> <layers:%d>\n", rank, layers); -> conversar com Victor sobre
if (dimensions == 2 && (sp_surface_tracking)) {
PetscPrintf(PETSC_COMM_WORLD,"\nSurface Processes Swarm (creating)\n");

sp_create_surface_swarm_2d();

// Surface Processes Swarm
if (dimensions == 2 && geoq_on && sp_surface_tracking && n_interfaces>0 && interfaces_from_ascii==1) {
PetscPrintf(PETSC_COMM_WORLD, "\nSP Swarm (creating)\n");
ierr = sp_create_surface_vec(); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "SP Swarm: done\n");
ierr = sp_interpolate_surface_particles_to_vec(); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"Surface Processes Swarm: done\n");
}

// PetscPrintf(PETSC_COMM_SELF,"********** <rank:%d> <particles_per_ele:%d>\n", rank, particles_per_ele); -> conversar com Victor sobre
// PetscPrintf(PETSC_COMM_SELF,"********** <rank:%d> <layers:%d>\n", rank, layers); -> conversar com Victor sobre

if (dimensions == 3) {
calc_drho(); // CHECK: does need this call here?
}
Expand Down Expand Up @@ -228,8 +229,9 @@ int main(int argc,char **args)
ierr = write_geoq_(tcont,binary_output);
ierr = write_tempo(tcont);

if (dimensions == 2 && sp_surface_tracking && geoq_on && n_interfaces>0 && interfaces_from_ascii==1) {
sp_write_surface_vec(tcont);
if (dimensions == 2 && sp_surface_tracking) {
ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN-1,"surface_%d", tcont); CHKERRQ(ierr);
ierr = sp_view_2d(dms_s, prefix); CHKERRQ(ierr);
}

VecCopy(Veloc_fut,Veloc);
Expand All @@ -242,14 +244,6 @@ int main(int argc,char **args)
PetscPrintf(PETSC_COMM_WORLD, "\n\n*** Using custom initial print_step = %d until %.3g Myr\n\n", print_step, initial_print_max_time);
}

if (dimensions == 2 && sp_surface_processes && PETSC_FALSE == set_sp_dt) {
sp_dt = 10.0 * dt_calor;
PetscPrintf(PETSC_COMM_WORLD, "\n\n*** Using default sp_dt\n\n");
}

sp_eval_time = sp_dt;
sp_last_eval_time = 0.0;

for (tempo = dt_calor,tcont=1;tempo<=timeMAX && tcont<=stepMAX;tempo+=dt_calor, tcont++){
if ((tempo > initial_print_max_time || fabs(tempo-initial_print_max_time) < 0.0001) && initial_print_step > 0) {
initial_print_step = 0;
Expand Down Expand Up @@ -277,17 +271,6 @@ int main(int argc,char **args)

ierr = veloc_total(dimensions); CHKERRQ(ierr);

if (dimensions == 2 && sp_surface_processes && geoq_on && n_interfaces>0 && interfaces_from_ascii==1 && (tempo > sp_eval_time || fabs(tempo-sp_eval_time) < 0.0001)) {
PetscPrintf(PETSC_COMM_WORLD,"\nEvaluating sp...\n");

ierr = rescalePrecipitation(tempo);

evaluate_surface_processes();

sp_eval_time += sp_dt;
sp_last_eval_time = tempo;
}

if (geoq_on){
if (RK4==1){
VecCopy(Veloc_fut,Veloc_weight);
Expand All @@ -302,6 +285,10 @@ int main(int argc,char **args)
// y a b x
VecAXPBY(Veloc_weight, fac, (1.0-fac),Veloc_fut); //y = a*x + b*y
ierr = moveSwarm(dimensions, dt_calor_sec/max_cont);

if (dimensions == 2 && sp_surface_tracking) {
ierr = sp_move_surface_swarm(dimensions, dt_calor_sec/max_cont); CHKERRQ(ierr);
}
}
}

Expand All @@ -312,8 +299,13 @@ int main(int argc,char **args)
}
}

if (dimensions == 2 && sp_surface_tracking && geoq_on && n_interfaces>0 && interfaces_from_ascii==1) {
ierr = sp_interpolate_surface_particles_to_vec(); CHKERRQ(ierr);
if (dimensions == 2 && sp_surface_tracking) {
ierr = sp_surface_swarm_interpolation(); CHKERRQ(ierr);

if (sp_surface_processes) {
ierr = sp_evaluate_surface_processes(dimensions, dt_calor_sec); CHKERRQ(ierr);
ierr = sp_update_surface_swarm_particles_properties(); CHKERRQ(ierr);
}
}

if (tcont%print_step==0){
Expand All @@ -334,8 +326,9 @@ int main(int argc,char **args)
}
}

if (dimensions == 2 && sp_surface_tracking && n_interfaces>0 && interfaces_from_ascii==1) {
ierr = sp_write_surface_vec(tcont); CHKERRQ(ierr);
if (dimensions == 2 && sp_surface_tracking) {
ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN-1,"surface_%d", tcont); CHKERRQ(ierr);
ierr = sp_view_2d(dms_s, prefix); CHKERRQ(ierr);
}
}
}
Expand All @@ -351,8 +344,7 @@ int main(int argc,char **args)

destroy_veloc();


if (sp_surface_processes && geoq_on && n_interfaces>0 && interfaces_from_ascii==1) {
if (dimensions == 2 && sp_surface_tracking) {
sp_destroy();
}

Expand Down Expand Up @@ -516,17 +508,3 @@ PetscErrorCode multi_veloc_change(Vec Veloc_fut,double tempo){
}
PetscFunctionReturn(0);
}


PetscErrorCode rescalePrecipitation(double tempo)
{
//PetscErrorCode ierr;
if (cont_var_climate<n_var_climate){
if (tempo>1.0E6*var_climate_time[cont_var_climate]){
prec_factor=var_climate_scale[cont_var_climate];
cont_var_climate++;
}
}

PetscFunctionReturn(0);
}
Loading

0 comments on commit 055b4b3

Please sign in to comment.