Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Shear heating #124

Merged
merged 4 commits into from
Feb 7, 2024
Merged

Shear heating #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
68fb80d
Implementation of shear heating
victorsacek Oct 16, 2023
a092732
Implementation of the shear heating for the 2D model
victorsacek Oct 16, 2023
9d36cad
correction in the term of shear heating: now in W/kg
victorsacek Oct 25, 2023
51e2852
Incorporation of Xi_shear to scale the shear heating
victorsacek Oct 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 17 additions & 5 deletions src/DMSwarm2mesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,9 @@ extern PetscScalar *inter_H;

extern PetscInt periodic_boundary;

extern PetscInt WITH_SHEAR_H;
extern PetscReal Xi_shear;

extern double visc_MIN;

extern PetscReal air_threshold_density;
Expand All @@ -43,6 +46,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 @@ -100,6 +105,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 @@ -149,33 +156,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
Open in desktop
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,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 @@ -265,6 +267,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 @@ -405,6 +409,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
Loading