diff --git a/src/DMSwarm2mesh.cpp b/src/DMSwarm2mesh.cpp index 59c2e924..56ac9252 100644 --- a/src/DMSwarm2mesh.cpp +++ b/src/DMSwarm2mesh.cpp @@ -3,6 +3,18 @@ #include PetscErrorCode mean_value_periodic_boundary_2d(DM da,Vec F,Vec local_F, PetscScalar **ff,int esc); +PetscScalar find_density(PetscScalar p, PetscScalar t, PetscInt idx); +PetscReal linear_interpolation(PetscReal rx, PetscReal rz,PetscScalar V0, PetscScalar V1, PetscScalar V2, PetscScalar V3); + +// void print_vec(PetscScalar *vec, PetscInt s_vec) +// { +// int i; +// for (i=0; i1) {printf("weird rx=%f , Swarm2Mesh\n",rx); exit(1);} if (rz<0 || rz>1) {printf("weird rz=%f , Swarm2Mesh\n",rz); exit(1);} + // Apply phase change density from ascii file + PetscReal rho_aux; + + tt_particle = linear_interpolation(rx,rz,tt[k][i],tt[k][i+1],tt[k+1][i],tt[k+1][i+1]); + pp_particle = linear_interpolation(rx,rz,pp_aux[k][i],pp_aux[k][i+1],pp_aux[k+1][i],pp_aux[k+1][i+1]); - + if (phase_change == 1 && phase_change_unit_number[layer_array[p]] >= 0) rho_aux = find_density(pp_particle, tt_particle, layer_array[p]); + else rho_aux = inter_rho[layer_array[p]] * (1.0 - alpha_exp_thermo * tt_particle); rfac = (1.0-rx)*(1.0-rz); - qq_rho [k][i] += rfac*inter_rho[layer_array[p]]; + qq_rho [k][i] += rfac*rho_aux; qq_H [k][i] += rfac*inter_H[layer_array[p]]; 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_rho [k][i+1] += rfac*rho_aux; qq_H [k][i+1] += rfac*inter_H[layer_array[p]]; 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_rho [k+1][i] += rfac*rho_aux; qq_H [k+1][i] += rfac*inter_H[layer_array[p]]; 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_rho [k+1][i+1] += rfac*rho_aux; qq_H [k+1][i+1] += rfac*inter_H[layer_array[p]]; qq_strain[k+1][i+1] += rfac*strain_fac[p]; qq_strain_rate[k+1][i+1] += rfac*strain_rate_fac[p]; @@ -319,6 +361,9 @@ 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_Temper,&tt);CHKERRQ(ierr); + ierr = DMDAVecRestoreArray(da_Thermal,local_P_aux,&pp_aux);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); @@ -406,7 +451,7 @@ PetscErrorCode Swarm2Mesh_2d(){ ierr = DMGlobalToLocalBegin(da_Thermal,Temper,INSERT_VALUES,local_Temper); ierr = DMGlobalToLocalEnd( da_Thermal,Temper,INSERT_VALUES,local_Temper); - ierr = DMDAVecGetArray(da_Thermal,local_Temper,&TT);CHKERRQ(ierr); + ierr = DMDAVecGetArray(da_Thermal,local_Temper,&tt);CHKERRQ(ierr); ierr = DMGlobalToLocalBegin(da_Thermal,geoq_rho,INSERT_VALUES,local_geoq_rho); ierr = DMGlobalToLocalEnd( da_Thermal,geoq_rho,INSERT_VALUES,local_geoq_rho); @@ -423,14 +468,14 @@ PetscErrorCode Swarm2Mesh_2d(){ for (k=sz; k phase_pressure[p_cum_size[file_index+1]-1]) p_value = phase_pressure[p_cum_size[file_index+1]-1]; + + if (t_value < phase_temperature[p_cum_size[file_index]]) t_value = phase_temperature[p_cum_size[file_index]]; + else if (t_value > phase_temperature[p_cum_size[file_index+1]-1]) t_value = phase_temperature[p_cum_size[file_index+1]-1]; + + p_index = p_cum_size[file_index] + round( + (p_cum_size[file_index+1]-p_cum_size[file_index]-1)* + (p_value-phase_pressure[p_cum_size[file_index]])/ + (phase_pressure[p_cum_size[file_index+1]-1]-phase_pressure[p_cum_size[file_index]])); + // Sanity check + if (p_index < p_cum_size[file_index]) p_index = p_cum_size[file_index]; + else if (p_index >= p_cum_size[file_index+1]) p_index = p_cum_size[file_index+1] - 1; + + t_index = t_cum_size[file_index] + round( + (t_cum_size[file_index+1]-t_cum_size[file_index]-1)* + (t_value-phase_temperature[t_cum_size[file_index]])/ + (phase_temperature[t_cum_size[file_index+1]-1]-phase_temperature[t_cum_size[file_index]])); + // Sanity check + if (t_index < t_cum_size[file_index]) t_index = t_cum_size[file_index]; + else if (t_index >= t_cum_size[file_index+1]) t_index = t_cum_size[file_index+1] - 1; + + d_index = (p_index-p_cum_size[file_index]) + d_cum_size[file_index] + ((t_index-t_cum_size[file_index])*(t_size[file_index+1])); + + // PetscPrintf(PETSC_COMM_WORLD, "p:%.5g, t:%.2f, layer_number:%d, file_idx:%d, idx_p:%d, idx_t:%d, idx_d:%d, density:%f\n", p_value, t_value, layer_number, file_index, p_index, t_index, d_index, phase_density[d_index]); + // PetscPrintf(PETSC_COMM_WORLD, "p_cum_size: %d\n", p_cum_size[file_index]); + // PetscPrintf(PETSC_COMM_WORLD, "p: %.5g -> p[%d]: %.5g\n", p_value, p_index, phase_pressure[p_index]); + // PetscPrintf(PETSC_COMM_WORLD, "t: %.5g -> t[%d]: %.5g\n", t_value, t_index, phase_temperature[t_index]); + // PetscPrintf(PETSC_COMM_WORLD, "d[%d]: %.5g\n", d_index, phase_density[d_index]); + PetscFunctionReturn(phase_density[d_index]); +} diff --git a/src/DMSwarm_move.cpp b/src/DMSwarm_move.cpp index 5865e1cb..d94085c1 100644 --- a/src/DMSwarm_move.cpp +++ b/src/DMSwarm_move.cpp @@ -5,6 +5,7 @@ #include #include #include "petscsys.h" +#include typedef struct { PetscScalar u; diff --git a/src/calc_drho.cpp b/src/calc_drho.cpp index 37dd51c9..3dcf9b07 100644 --- a/src/calc_drho.cpp +++ b/src/calc_drho.cpp @@ -13,13 +13,17 @@ PetscErrorCode calc_drho() { PetscErrorCode ierr=0; - VecCopy(Temper, dRho); CHKERRQ(ierr); + VecCopy(geoq_rho,dRho); - PetscScalar a; + VecScale(dRho,-gravity); - a=alpha_exp_thermo*gravity*RHOM; CHKERRQ(ierr);//check: replace RHOM by geoq_rho + // VecCopy(Temper, dRho); CHKERRQ(ierr); - VecAXPBY(dRho,-gravity,a,geoq_rho); CHKERRQ(ierr); + // PetscScalar a; + + // a=alpha_exp_thermo*gravity*RHOM; CHKERRQ(ierr);//check: replace RHOM by geoq_rho + + // VecAXPBY(dRho,-gravity,a,geoq_rho); CHKERRQ(ierr); return ierr; diff --git a/src/header.h b/src/header.h index 6b16d6bc..fea3082d 100644 --- a/src/header.h +++ b/src/header.h @@ -49,6 +49,7 @@ PetscInt checkered = 0; // 1=True, 0=False PetscInt initial_dynamic_range = 0; // 1=True, 0=False PetscInt periodic_boundary = 0; // 1=True, 0=False PetscInt high_kappa_in_asthenosphere = 0; // 1=True, 0=False +PetscInt phase_change = 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 @@ -98,6 +99,19 @@ int bcT_left; int bcT_right; // End of parameter file variables +// Phase change variables +PetscInt num_phase_change_files = 0; +PetscInt *phase_change_unit_number; +PetscInt *p_size, *p_size_0; +PetscInt *p_cum_size, *p_cum_size_0; +PetscInt *t_size, *t_size_0; +PetscInt *t_cum_size, *t_cum_size_0; +PetscInt *d_size, *d_size_0; +PetscInt *d_cum_size, *d_cum_size_0; +PetscScalar *phase_pressure, *phase_pressure_0; +PetscScalar *phase_temperature, *phase_temperature_0; +PetscScalar *phase_density, *phase_density_0; + PetscInt Px = PETSC_DECIDE; PetscInt Py = PETSC_DECIDE; PetscInt Pz = PETSC_DECIDE; @@ -386,4 +400,4 @@ PetscReal adiabatic_scaled; PetscReal strain_rate0_scaled; PetscReal pressure0_scaled; -PetscReal air_threshold_density; +PetscReal air_threshold_density; \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 95ce7e98..d508975a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -328,6 +328,20 @@ int main(int argc,char **args) PetscTime(&Tempo2); + if (phase_change==1) + { + PetscFree(phase_pressure); + PetscFree(phase_temperature); + PetscFree(phase_density); + PetscFree(phase_change_unit_number); + PetscFree(p_size); + PetscFree(p_cum_size); + PetscFree(t_size); + PetscFree(t_cum_size); + PetscFree(d_size); + PetscFree(d_cum_size); + } + ierr = PetscFinalize(); return 0; diff --git a/src/reader.cpp b/src/reader.cpp index 2d619626..6de85ce9 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -1,10 +1,33 @@ #include +#include // Prototypes int check_a_b(char tkn_w[], char tkn_v[], const char str_a[], const char str_b[]); PetscBool check_a_b_bool(char tkn_w[], char tkn_v[], const char str_a[], const char str_b[]); //void ErrorInterfaces(int rank, const char fname[], int flag); void ErrorInterfaces(); +PetscInt read_phase_change_file(char *fname, PetscInt file_index); +PetscBool is_positive_number(char tkn); + +void print_vec(PetscScalar *vec, PetscInt s_vec) +{ + int i; + for (i=0; i\n", fname); + if (file_index == 0) + { + PetscCalloc1(file_index+2,&p_size_0); + PetscCalloc1(file_index+2,&p_cum_size_0); + PetscCalloc1(file_index+2,&t_size_0); + PetscCalloc1(file_index+2,&t_cum_size_0); + PetscCalloc1(file_index+2,&d_size_0); + PetscCalloc1(file_index+2,&d_cum_size_0); + } + else + { + PetscRealloc((file_index+2)*sizeof(PetscInt),&p_size_0); + PetscRealloc((file_index+2)*sizeof(PetscInt),&p_cum_size_0); + PetscRealloc((file_index+2)*sizeof(PetscInt),&t_size_0); + PetscRealloc((file_index+2)*sizeof(PetscInt),&t_cum_size_0); + PetscRealloc((file_index+2)*sizeof(PetscInt),&d_size_0); + PetscRealloc((file_index+2)*sizeof(PetscInt),&d_cum_size_0); + } + + while (!feof(file)) + { + nread = fgets(line, size, file); + if ((nread == NULL) || (strlen(line) == 0) || (line[0] == '\n') || (line[0] == '#')) continue; + tkn_0 = strtok(line, " \t\n"); + if (tkn_0[0] == '#') continue; + if (strcmp(tkn_0, "unit") == 0) + { + for (i=0; i on file <%s>.\n", tkn_0, fname); + exit(1); + } + if (tkn_1 == NULL) break; + phase_change_unit_number[atoi(tkn_1)] = file_index; + } + } + else if (strcmp(tkn_0, "pressure") == 0) + { + for (i=0; i<3; i+=1) + { + tkn_1 = strtok(NULL, " \t\n"); + if (tkn_1 == NULL) + { + PetscPrintf(PETSC_COMM_WORLD, "Error. Insufficient values for <%s> on file <%s>.\n", tkn_0, fname); + exit(1); + } + if (i<2) {aux[i] = atof(tkn_1);} + else + { + p_size_0[file_index+1] = atoi(tkn_1); + p_cum_size_0[file_index+1] = p_size_0[file_index+1] + p_cum_size_0[file_index]; + } + } + + if (file_index == 0) PetscCalloc1(p_cum_size_0[file_index+1],&phase_pressure_0); + else PetscRealloc(p_cum_size_0[file_index+1]*sizeof(PetscScalar), &phase_pressure_0); + + for (i=p_cum_size_0[file_index]; i on file <%s>.\n", tkn_0, fname); + exit(1); + } + if (i<2) {aux[i] = atof(tkn_1);} + else + { + t_size_0[file_index+1] = atoi(tkn_1); + t_cum_size_0[file_index+1] = t_size_0[file_index+1] + t_cum_size_0[file_index]; + } + } + + if (file_index == 0) PetscCalloc1(t_cum_size_0[file_index+1],&phase_temperature_0); + else PetscRealloc(t_cum_size_0[file_index+1]*sizeof(PetscScalar), &phase_temperature_0); + + for (i=t_cum_size_0[file_index]; i on file <%s>.\n", fname); + exit(1); + } + ix = i + corr_i + (j * corr_j); + phase_density_0[ix] = atof(tkn_0); + tkn_0 = strtok(NULL, " \t\n"); + } + if (j>=t_size_0[file_index+1]) + { + PetscPrintf(PETSC_COMM_WORLD, "Error. Too many values for on file <%s>.\n", fname); + exit(1); + } + j += 1; + } + else + { + PetscPrintf(PETSC_COMM_WORLD, "Error. Unrecognized keyword/value <%s> on file <%s>.\n", tkn_0, fname); + exit(1); + } + } + fclose(file); + + if (j on file <%s>.\n", fname); + exit(1); + } + PetscFunctionReturn(0); +} + +PetscBool is_positive_number(char tkn) +{ + if (tkn=='0'||tkn=='1'||tkn=='2'||tkn=='3'||tkn=='4'||tkn=='5'||tkn=='6'||tkn=='7'||tkn=='8'||tkn=='9') PetscFunctionReturn(PETSC_TRUE); + else PetscFunctionReturn(PETSC_FALSE); +} +