Skip to content

Commit

Permalink
Physical parameters - thermal conductivity variable by lithology (#121)
Browse files Browse the repository at this point in the history 
* first commit

* Implementation of thermal conductivity function of the lithology: complete

* correction in the code

* The option viscosity_per_element = variable is no longer available.

* thermal diffusivity constant per element

* export thermal diffusivity

* flag export_thermal_diffusivity created. Default = False

* cleaning of the code

* example of numerical model with two layers with different conductivity.

* Modification of the interfaces.txt file to the original values

* add test for variable conductivity

---------

Co-authored-by: Jamison Assunção <[email protected]>
Co-authored-by: Rafael Silva <[email protected]>
  • Loading branch information
@victorsacek @jamisonassuncao @rafaelmds
3 people committed Nov 8, 2023
1 parent 1d4b12e commit 19cc610
Show file tree
Hide file tree
Showing 35 changed files with 3,790 additions and 129 deletions.
763 changes: 763 additions & 0 deletions examples/two_layers_with_variable_conductivity/create_input_and_plot.ipynb
View file

Large diffs are not rendered by default.

31 changes: 29 additions & 2 deletions src/DM1_2d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ extern Vec local_FT;
extern Vec local_Temper;

extern Vec geoq_H,local_geoq_H;
extern Vec geoq_kappa, local_geoq_kappa;

extern Vec local_TC;

Expand Down Expand Up @@ -176,7 +177,17 @@ PetscErrorCode AssembleA_Thermal_2d(Mat A,DM thermal_da,PetscReal *TKe,PetscReal
ierr = DMGlobalToLocalEnd(thermal_da,Temper,INSERT_VALUES,local_Temper);

ierr = DMDAVecGetArray(thermal_da,local_Temper,&tt);CHKERRQ(ierr);

////////

PetscScalar **qq_kappa;

ierr = VecZeroEntries(local_geoq_kappa);CHKERRQ(ierr);

ierr = DMGlobalToLocalBegin(thermal_da,geoq_kappa,INSERT_VALUES,local_geoq_kappa);
ierr = DMGlobalToLocalEnd(thermal_da,geoq_kappa,INSERT_VALUES,local_geoq_kappa);

ierr = DMDAVecGetArray(thermal_da,local_geoq_kappa,&qq_kappa);CHKERRQ(ierr);


PetscInt M,P;
Expand Down Expand Up @@ -240,7 +251,9 @@ PetscErrorCode AssembleA_Thermal_2d(Mat A,DM thermal_da,PetscReal *TKe,PetscReal
v_vec_aux_ele[i*2+1] = VV[ind[i].j][ind[i].i].w;
}

kappa_eff = kappa;
kappa_eff = qq_kappa[ek][ei];

//kappa_eff = kappa;
if (high_kappa_in_asthenosphere==1){
tt_ele[0] = tt[ek][ei];
tt_ele[1] = tt[ek][ei+1];
Expand Down Expand Up @@ -287,6 +300,7 @@ PetscErrorCode AssembleA_Thermal_2d(Mat A,DM thermal_da,PetscReal *TKe,PetscReal
ierr = DMDAVecRestoreArray(veloc_da,local_V,&VV);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(thermal_da,local_TC,&TTC);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(thermal_da,local_Temper,&tt);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(thermal_da,local_geoq_kappa,&qq_kappa);CHKERRQ(ierr);


PetscFunctionReturn(0);
Expand Down Expand Up @@ -335,6 +349,16 @@ PetscErrorCode AssembleF_Thermal_2d(Vec F,DM thermal_da,PetscReal *TKe,PetscReal
////////


PetscScalar **qq_kappa;

ierr = VecZeroEntries(local_geoq_kappa);CHKERRQ(ierr);

ierr = DMGlobalToLocalBegin(thermal_da,geoq_kappa,INSERT_VALUES,local_geoq_kappa);
ierr = DMGlobalToLocalEnd(thermal_da,geoq_kappa,INSERT_VALUES,local_geoq_kappa);

ierr = DMDAVecGetArray(thermal_da,local_geoq_kappa,&qq_kappa);CHKERRQ(ierr);



/* get acces to the vector */

Expand Down Expand Up @@ -388,7 +412,9 @@ PetscErrorCode AssembleF_Thermal_2d(Vec F,DM thermal_da,PetscReal *TKe,PetscReal
v_vec_aux_ele[i*2+1] = VV[ind[i].j][ind[i].i].w;
}

kappa_eff = kappa;
kappa_eff = qq_kappa[ek][ei];

//kappa_eff = kappa;
if (high_kappa_in_asthenosphere==1){
tt_ele[0] = tt[ek][ei];
tt_ele[1] = tt[ek][ei+1];
Expand Down Expand Up @@ -486,6 +512,7 @@ PetscErrorCode AssembleF_Thermal_2d(Vec F,DM thermal_da,PetscReal *TKe,PetscReal
ierr = DMDAVecRestoreArray(thermal_da,local_TC,&TTC);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(thermal_da,local_Temper,&tt);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(thermal_da,local_geoq_H,&HH);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(thermal_da,local_geoq_kappa,&qq_kappa);CHKERRQ(ierr);

if (periodic_boundary==1){

Expand Down
25 changes: 22 additions & 3 deletions src/DMSwarm2mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ extern Vec geoq_cont,local_geoq_cont;
extern Vec geoq_H,local_geoq_H;
extern Vec geoq_strain, local_geoq_strain;
extern Vec geoq_strain_rate, local_geoq_strain_rate;
extern Vec geoq_kappa, local_geoq_kappa;

extern Vec Temper,local_Temper;

Expand All @@ -32,6 +33,7 @@ extern PetscInt sticky_blanket_air;

extern PetscScalar *inter_rho;
extern PetscScalar *inter_H;
extern PetscScalar *conductivity;

extern PetscInt periodic_boundary;

Expand All @@ -46,7 +48,7 @@ extern PetscReal epsilon_x;
PetscErrorCode Swarm2Mesh_2d(){

PetscErrorCode ierr;
PetscScalar **qq,**qq_cont,**qq_rho,**TT,**qq_H,**qq_strain;
PetscScalar **qq,**qq_cont,**qq_rho,**TT,**qq_H,**qq_strain,**qq_kappa;
PetscScalar **qq_strain_rate;

ierr = VecSet(geoq,0.0);CHKERRQ(ierr);
Expand All @@ -55,12 +57,14 @@ PetscErrorCode Swarm2Mesh_2d(){
ierr = VecSet(geoq_H,0.0);CHKERRQ(ierr);
ierr = VecSet(geoq_strain,0.0);CHKERRQ(ierr);
ierr = VecSet(geoq_strain_rate,0.0);CHKERRQ(ierr);


ierr = VecZeroEntries(local_geoq);CHKERRQ(ierr);
ierr = VecZeroEntries(local_geoq_rho);CHKERRQ(ierr);
ierr = VecZeroEntries(local_geoq_H);CHKERRQ(ierr);
ierr = VecZeroEntries(local_geoq_strain);CHKERRQ(ierr);
ierr = VecZeroEntries(local_geoq_strain_rate);CHKERRQ(ierr);


ierr = DMGlobalToLocalBegin(da_Thermal,geoq,INSERT_VALUES,local_geoq);
ierr = DMGlobalToLocalEnd( da_Thermal,geoq,INSERT_VALUES,local_geoq);
Expand All @@ -77,7 +81,6 @@ PetscErrorCode Swarm2Mesh_2d(){
ierr = DMGlobalToLocalBegin(da_Thermal,geoq_strain_rate,INSERT_VALUES,local_geoq_strain_rate);
ierr = DMGlobalToLocalEnd( da_Thermal,geoq_strain_rate,INSERT_VALUES,local_geoq_strain_rate);


ierr = DMDAVecGetArray(da_Thermal,local_geoq,&qq);CHKERRQ(ierr);
ierr = DMDAVecGetArray(da_Thermal,local_geoq_rho,&qq_rho);CHKERRQ(ierr);
ierr = DMDAVecGetArray(da_Thermal,local_geoq_H,&qq_H);CHKERRQ(ierr);
Expand Down Expand Up @@ -183,6 +186,9 @@ PetscErrorCode Swarm2Mesh_2d(){

if (visc_const_per_element==0){

PetscPrintf(PETSC_COMM_WORLD,"\n\nThe option \nviscosity_per_element = variable \nis no longer available.\n");
exit(-2);

if (visc_harmonic_mean==1){
for (p=0; p<nlocal; p++) {
PetscReal cx,cz;
Expand Down Expand Up @@ -319,6 +325,7 @@ PetscErrorCode Swarm2Mesh_2d(){
ierr = DMLocalToGlobalBegin(da_Thermal,local_geoq_cont,ADD_VALUES,geoq_cont);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(da_Thermal,local_geoq_cont,ADD_VALUES,geoq_cont);CHKERRQ(ierr);


if (periodic_boundary==1){
ierr = mean_value_periodic_boundary_2d(da_Thermal,geoq_rho,local_geoq_rho,qq_rho,1);
ierr = mean_value_periodic_boundary_2d(da_Thermal,geoq_H,local_geoq_H,qq_H,1);
Expand All @@ -340,11 +347,15 @@ PetscErrorCode Swarm2Mesh_2d(){
if (visc_const_per_element==1){
ierr = VecSet(geoq,0.0);CHKERRQ(ierr);
ierr = VecSet(geoq_cont,0.0);CHKERRQ(ierr);
ierr = VecSet(geoq_kappa,0.0);CHKERRQ(ierr);

ierr = DMGlobalToLocalBegin(da_Thermal,geoq,INSERT_VALUES,local_geoq);
ierr = DMGlobalToLocalEnd( da_Thermal,geoq,INSERT_VALUES,local_geoq);
ierr = DMDAVecGetArray(da_Thermal,local_geoq,&qq);CHKERRQ(ierr);

ierr = DMGlobalToLocalBegin(da_Thermal,geoq_kappa,INSERT_VALUES,local_geoq_kappa);
ierr = DMGlobalToLocalEnd( da_Thermal,geoq_kappa,INSERT_VALUES,local_geoq_kappa);
ierr = DMDAVecGetArray(da_Thermal,local_geoq_kappa,&qq_kappa);CHKERRQ(ierr);

ierr = DMGlobalToLocalBegin(da_Thermal,geoq_cont,INSERT_VALUES,local_geoq_cont);
ierr = DMGlobalToLocalEnd( da_Thermal,geoq_cont,INSERT_VALUES,local_geoq_cont);
Expand All @@ -361,6 +372,7 @@ PetscErrorCode Swarm2Mesh_2d(){
k = (int)((cz+depth)/dz_const);

qq[k][i] += 1.0/geoq_fac[p];
qq_kappa[k][i] += 1.0/conductivity[layer_array[p]];
qq_cont[k][i] += 1.0;
}
}
Expand All @@ -376,6 +388,7 @@ PetscErrorCode Swarm2Mesh_2d(){
k = (int)((cz+depth)/dz_const);

qq[k][i] += geoq_fac[p];
qq_kappa[k][i] += conductivity[layer_array[p]];
qq_cont[k][i] += 1.0;
}
}
Expand All @@ -389,8 +402,14 @@ PetscErrorCode Swarm2Mesh_2d(){
ierr = DMLocalToGlobalBegin(da_Thermal,local_geoq_cont,ADD_VALUES,geoq_cont);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(da_Thermal,local_geoq_cont,ADD_VALUES,geoq_cont);CHKERRQ(ierr);

ierr = DMDAVecRestoreArray(da_Thermal,local_geoq_kappa,&qq_kappa);CHKERRQ(ierr);
ierr = DMLocalToGlobalBegin(da_Thermal,local_geoq_kappa,ADD_VALUES,geoq_kappa);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(da_Thermal,local_geoq_kappa,ADD_VALUES,geoq_kappa);CHKERRQ(ierr);

VecPointwiseDivide(geoq,geoq,geoq_cont);
VecPointwiseDivide(geoq_kappa,geoq_kappa,geoq_cont);
if (visc_harmonic_mean==1) VecReciprocal(geoq); //<--- harmonic
if (visc_harmonic_mean==1) VecReciprocal(geoq_kappa); //<--- harmonic

}

Expand Down Expand Up @@ -470,7 +489,7 @@ PetscErrorCode Swarm2Mesh_2d(){

}

PetscErrorCode Swarm2Mesh_3d(){
PetscErrorCode Swarm2Mesh_3d(){ ///!!! geoq_kappa for 3D models have to be implemented!!!

PetscErrorCode ierr;
PetscScalar ***qq,***qq_cont,***qq_rho,***TT,***qq_H,***qq_strain;
Expand Down
10 changes: 10 additions & 0 deletions src/DMT.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ extern Vec geoq_rho,local_geoq_rho;
extern Vec geoq_H,local_geoq_H;
extern Vec geoq_strain,local_geoq_strain;
extern Vec geoq_strain_rate,local_geoq_strain_rate;
extern Vec geoq_kappa, local_geoq_kappa;

extern Vec geoq_cont,local_geoq_cont;

Expand Down Expand Up @@ -101,6 +102,8 @@ extern double kappa;
extern double RHOM;
extern double c_heat_capacity;

extern PetscBool export_kappa;

PetscErrorCode create_thermal(int dimensions, PetscInt mx, PetscInt my, PetscInt mz, PetscInt Px, PetscInt Py, PetscInt Pz)
{

Expand Down Expand Up @@ -201,6 +204,7 @@ PetscErrorCode create_thermal(int dimensions, PetscInt mx, PetscInt my, PetscInt
ierr = DMCreateGlobalVector(da_Thermal,&geoq_strain);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da_Thermal,&geoq_strain_rate);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da_Thermal,&geoq_cont);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(da_Thermal,&geoq_kappa);CHKERRQ(ierr);

ierr = DMCreateGlobalVector(da_Thermal,&dRho);CHKERRQ(ierr);

Expand All @@ -222,6 +226,7 @@ PetscErrorCode create_thermal(int dimensions, PetscInt mx, PetscInt my, PetscInt
ierr = DMCreateLocalVector(da_Thermal,&local_geoq_strain);
ierr = DMCreateLocalVector(da_Thermal,&local_geoq_strain_rate);
ierr = DMCreateLocalVector(da_Thermal,&local_geoq_cont);
ierr = DMCreateLocalVector(da_Thermal,&local_geoq_kappa);

ierr = DMCreateLocalVector(da_Thermal,&local_dRho);

Expand Down Expand Up @@ -530,5 +535,10 @@ PetscErrorCode write_geoq_(int cont, PetscInt binary_out)
sprintf(variable_name,"strain_rate");
write_all_(cont,geoq_strain_rate,variable_name,binary_out);

if (export_kappa==PETSC_TRUE){
sprintf(variable_name,"thermal_diffusivity");
write_all_(cont,geoq_kappa,variable_name,binary_out);
}

PetscFunctionReturn(0);
}
5 changes: 5 additions & 0 deletions src/header.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,6 +142,7 @@ PetscScalar *cohesion_max;
PetscScalar *friction_angle_min;
PetscScalar *friction_angle_max;
PetscBool weakening_from_interfaces_file = PETSC_FALSE;
PetscScalar *conductivity;

int tcont=0;

Expand Down Expand Up @@ -234,6 +235,8 @@ Vec local_geoq_cont;
Vec geoq_strain;
Vec local_geoq_strain;

Vec geoq_kappa;
Vec local_geoq_kappa;

Vec geoq_strain_rate;
Vec local_geoq_strain_rate;
Expand Down Expand Up @@ -405,3 +408,5 @@ PetscReal strain_rate0_scaled;
PetscReal pressure0_scaled;

PetscReal air_threshold_density;

PetscBool export_kappa = PETSC_FALSE;
Loading

0 comments on commit 19cc610

Please sign in to comment.