Skip to content

Commit

Permalink
Shear heating (#124)
Browse files Browse the repository at this point in the history 
* Implementation of shear heating

* Implementation of the shear heating for the 2D model

* correction in the term of shear heating: now in W/kg

* Incorporation of Xi_shear to scale the shear heating
  • Loading branch information
@victorsacek
victorsacek authored Feb 7, 2024
1 parent 19cc610 commit 6089f97
Show file tree
Hide file tree
Showing 3 changed files with 25 additions and 5 deletions.
22 changes: 17 additions & 5 deletions src/DMSwarm2mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@ extern PetscScalar *conductivity;

extern PetscInt periodic_boundary;

extern PetscInt WITH_SHEAR_H;
extern PetscReal Xi_shear;

extern double visc_MIN;

extern PetscReal air_threshold_density;
Expand All @@ -45,6 +48,8 @@ extern PetscReal rho0_scaled;

extern PetscReal epsilon_x;

//Shear heating only in 2D!!! Must be implemented for 3D!!!

PetscErrorCode Swarm2Mesh_2d(){

PetscErrorCode ierr;
Expand Down Expand Up @@ -103,6 +108,8 @@ PetscErrorCode Swarm2Mesh_2d(){
PetscReal *strain_rate_fac;
PetscInt *layer_array;

PetscReal inter_H_aux;

ierr = DMSwarmGetLocalSize(dms,&nlocal);CHKERRQ(ierr);

ierr = DMSwarmGetField(dms,DMSwarmPICField_coor,&bs,NULL,(void**)&array);CHKERRQ(ierr);
Expand Down Expand Up @@ -152,33 +159,38 @@ PetscErrorCode Swarm2Mesh_2d(){
if (rx<0 || rx>1) {printf("weird rx=%f , Swarm2Mesh\n",rx); exit(1);}
if (rz<0 || rz>1) {printf("weird rz=%f , Swarm2Mesh\n",rz); exit(1);}


if (WITH_SHEAR_H == 1){
inter_H_aux = inter_H[layer_array[p]] + Xi_shear*4*geoq_fac[p]*strain_rate_fac[p]*strain_rate_fac[p]/inter_rho[layer_array[p]];
}
else {
inter_H_aux = inter_H[layer_array[p]];
}


rfac = (1.0-rx)*(1.0-rz);
qq_rho [k][i] += rfac*inter_rho[layer_array[p]];
qq_H [k][i] += rfac*inter_H[layer_array[p]];
qq_H [k][i] += rfac*inter_H_aux;
qq_strain[k][i] += rfac*strain_fac[p];
qq_strain_rate[k][i] += rfac*strain_rate_fac[p];
qq_cont [k][i] += rfac;

rfac = (rx)*(1.0-rz);
qq_rho [k][i+1] += rfac*inter_rho[layer_array[p]];
qq_H [k][i+1] += rfac*inter_H[layer_array[p]];
qq_H [k][i+1] += rfac*inter_H_aux;
qq_strain[k][i+1] += rfac*strain_fac[p];
qq_strain_rate[k][i+1] += rfac*strain_rate_fac[p];
qq_cont [k][i+1] += rfac;

rfac = (1.0-rx)*(rz);
qq_rho [k+1][i] += rfac*inter_rho[layer_array[p]];
qq_H [k+1][i] += rfac*inter_H[layer_array[p]];
qq_H [k+1][i] += rfac*inter_H_aux;
qq_strain[k+1][i] += rfac*strain_fac[p];
qq_strain_rate[k+1][i] += rfac*strain_rate_fac[p];
qq_cont [k+1][i] += rfac;

rfac = (rx)*(rz);
qq_rho [k+1][i+1] += rfac*inter_rho[layer_array[p]];
qq_H [k+1][i+1] += rfac*inter_H[layer_array[p]];
qq_H [k+1][i+1] += rfac*inter_H_aux;
qq_strain[k+1][i+1] += rfac*strain_fac[p];
qq_strain_rate[k+1][i+1] += rfac*strain_rate_fac[p];
qq_cont [k+1][i+1] += rfac;
Expand Down
2 changes: 2 additions & 0 deletions src/header.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ PetscScalar sp_d_c = 0.0;
// Parameter file boolean variables
PetscInt WITH_NON_LINEAR = 0; // 1=True, 0=False
PetscInt WITH_ADIABATIC_H = 0; // 1=True, 0=False
PetscInt WITH_SHEAR_H = 0; // 1=True, 0=False
PetscReal Xi_shear = 1; //scale factor
PetscInt WITH_RADIOGENIC_H = 0; // 1=True, 0=False
PetscInt PLASTICITY = 1; // 1=True, 0=False
PetscInt direct_solver = 1; // 1=direct, 0=iterative
Expand Down
6 changes: 6 additions & 0 deletions src/reader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,8 @@ extern int bcT_front;
extern int bcT_back;
extern PetscInt WITH_NON_LINEAR;
extern PetscInt WITH_ADIABATIC_H;
extern PetscInt WITH_SHEAR_H;
extern PetscReal Xi_shear;
extern PetscInt WITH_RADIOGENIC_H;
extern PetscInt PLASTICITY;
extern PetscReal denok_min;
Expand Down Expand Up @@ -268,6 +270,8 @@ PetscErrorCode reader(int rank, const char fName[]){
else if (strcmp(tkn_w, "geoq") == 0) {geoq_on = check_a_b(tkn_w, tkn_v, "on", "off");}
else if (strcmp(tkn_w, "non_linear_method") == 0) {WITH_NON_LINEAR = check_a_b(tkn_w, tkn_v, "on", "off");}
else if (strcmp(tkn_w, "adiabatic_component") == 0) {WITH_ADIABATIC_H = check_a_b(tkn_w, tkn_v, "on", "off");}
else if (strcmp(tkn_w, "shear_heating") == 0) {WITH_SHEAR_H = check_a_b(tkn_w, tkn_v, "on", "off");}
else if (strcmp(tkn_w, "Xi_shear") == 0) {Xi_shear = atof(tkn_v);}
else if (strcmp(tkn_w, "radiogenic_component") == 0) {WITH_RADIOGENIC_H = check_a_b(tkn_w, tkn_v, "on", "off");}
else if (strcmp(tkn_w, "plasticity") == 0) {PLASTICITY = check_a_b(tkn_w, tkn_v, "on", "off");}
else if (strcmp(tkn_w, "top_normal_velocity") == 0) {bcv_top_normal = check_a_b(tkn_w, tkn_v, "fixed", "free");}
Expand Down Expand Up @@ -409,6 +413,8 @@ PetscErrorCode reader(int rank, const char fName[]){
MPI_Bcast(&c_heat_capacity,1,MPI_DOUBLE,0,PETSC_COMM_WORLD);
MPI_Bcast(&WITH_NON_LINEAR,1,MPI_INT,0,PETSC_COMM_WORLD);
MPI_Bcast(&WITH_ADIABATIC_H,1,MPI_INT,0,PETSC_COMM_WORLD);
MPI_Bcast(&WITH_SHEAR_H,1,MPI_INT,0,PETSC_COMM_WORLD);
MPI_Bcast(&Xi_shear,1,MPIU_REAL,0,PETSC_COMM_WORLD);
MPI_Bcast(&WITH_RADIOGENIC_H,1,MPI_INT,0,PETSC_COMM_WORLD);
MPI_Bcast(&PLASTICITY,1,MPI_INT,0,PETSC_COMM_WORLD);
MPI_Bcast(&bcv_top_normal,1,MPI_INT,0,PETSC_COMM_WORLD);
Expand Down

0 comments on commit 6089f97

Please sign in to comment.