add new surface processes implementation (#125)

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`.

Author: rafaelmds

src/header.h

@@ -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;
@@ -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
@@ -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!)
@@ -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
@@ -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)
@@ -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

src/main.cpp

@@ -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)
 {
@@ -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
@@ -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?
 	}
@@ -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);
@@ -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;
@@ -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);
@@ -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);
+					}
 				}
 			}
 
@@ -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){
@@ -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);
 				}
 			}
 		}
@@ -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();
 	}
 
@@ -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);
-}