Skip to content

Commit

Permalink
diode tests
Browse files Browse the repository at this point in the history
wortiz committed Mar 5, 2024
1 parent b408fc3 commit c438c17
Showing 3 changed files with 129 additions and 14 deletions.
38 changes: 30 additions & 8 deletions src/mm_fill_porous.c
Original file line number Diff line number Diff line change
@@ -6469,8 +6469,18 @@ double load_cap_pres(int ipore, int ilnode, int ignode, double saturation)
m_inv = -1.0 / mexp;

/* Normalized saturation - should be between 0 to 1 */
sat_norm = (saturation - sat_min) / (sat_max - sat_min);
d_sat_norm_d_saturation = 1.0 / (sat_max - sat_min);
if (saturation < (1e-9 * (sat_max - sat_min) + sat_min)) {
sat_norm = 1e-9;
d_sat_norm_d_saturation = 0.0;
} else if (saturation > sat_max) {
sat_norm = 1.0;
d_sat_norm_d_saturation = 0.0;
} else {
sat_norm = (saturation - sat_min) / (sat_max - sat_min);
d_sat_norm_d_saturation = 1.0 / (sat_max - sat_min);
}
// sat_norm = (saturation - sat_min) / (sat_max - sat_min);
// d_sat_norm_d_saturation = 1.0 / (sat_max - sat_min);

brack_in = pow(sat_norm, m_inv) - 1.0;
d_brack_in_d_sat_norm = m_inv * pow(sat_norm, (m_inv - 1.0));
@@ -6490,18 +6500,30 @@ double load_cap_pres(int ipore, int ilnode, int ignode, double saturation)

switch (ipore) {
case 0:
mp->d_cap_pres[SHELL_SAT_1] = n_inv * con_c * pow(brack_in, (n_inv - 1.0)) *
d_brack_in_d_sat_norm * d_sat_norm_d_saturation;
if (brack_in > 0.0) {
mp->d_cap_pres[SHELL_SAT_1] = n_inv * con_c * pow(brack_in, (n_inv - 1.0)) *
d_brack_in_d_sat_norm * d_sat_norm_d_saturation;
} else {
mp->d_cap_pres[SHELL_SAT_1] = 0.0;
}
break;

case 1:
mp->d_cap_pres[SHELL_SAT_2] = n_inv * con_c * pow(brack_in, (n_inv - 1.0)) *
d_brack_in_d_sat_norm * d_sat_norm_d_saturation;
if (brack_in > 0.0) {
mp->d_cap_pres[SHELL_SAT_2] = n_inv * con_c * pow(brack_in, (n_inv - 1.0)) *
d_brack_in_d_sat_norm * d_sat_norm_d_saturation;
} else {
mp->d_cap_pres[SHELL_SAT_2] = 0.0;
}
break;

case 2:
mp->d_cap_pres[SHELL_SAT_3] = n_inv * con_c * pow(brack_in, (n_inv - 1.0)) *
d_brack_in_d_sat_norm * d_sat_norm_d_saturation;
if (brack_in > 0.0) {
mp->d_cap_pres[SHELL_SAT_3] = n_inv * con_c * pow(brack_in, (n_inv - 1.0)) *
d_brack_in_d_sat_norm * d_sat_norm_d_saturation;
} else {
mp->d_cap_pres[SHELL_SAT_3] = 0.0;
}
break;

default:
10 changes: 10 additions & 0 deletions src/mm_fill_shell.c
Original file line number Diff line number Diff line change
@@ -13265,6 +13265,7 @@ int assemble_porous_shell_saturation(dbl tt, // Time integration form

// Group saturation values at nodes and Gauss point
dbl sat_nodes[MAX_POR_SHELL][MDE] = {{0.0}};
dbl sat_nodes_old[MAX_POR_SHELL][MDE] = {{0.0}};

Check warning on line 13268 in src/mm_fill_shell.c

GitHub Actions / build release unittest

variable ‘sat_nodes_old’ set but not used [-Wunused-but-set-variable]
dbl sat_dot_nodes[MAX_POR_SHELL][MDE] = {{0.0}};
dbl sat_gauss[MAX_POR_SHELL] = {0.0};
dbl grad_sat_gauss[MAX_POR_SHELL][DIM] = {{0.0}};
@@ -13277,14 +13278,17 @@ int assemble_porous_shell_saturation(dbl tt, // Time integration form
switch (ipore) {
case 0:
sat_nodes[ipore][j] = *esp->sh_sat_1[j];
sat_nodes_old[ipore][j] = *esp_old->sh_sat_1[j];
sat_dot_nodes[ipore][j] = *esp_dot->sh_sat_1[j];
break;
case 1:
sat_nodes[ipore][j] = *esp->sh_sat_2[j];
sat_nodes_old[ipore][j] = *esp_old->sh_sat_2[j];
sat_dot_nodes[ipore][j] = *esp_dot->sh_sat_2[j];
break;
case 2:
sat_nodes[ipore][j] = *esp->sh_sat_3[j];
sat_nodes_old[ipore][j] = *esp_old->sh_sat_3[j];
sat_dot_nodes[ipore][j] = *esp_dot->sh_sat_3[j];
break;
}
@@ -13379,8 +13383,14 @@ int assemble_porous_shell_saturation(dbl tt, // Time integration form
if (pd->v[pg->imtrx][var]) {
/* Old method */
for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
// #define LAG_CAP_PRESSURE
#ifdef LAG_CAP_PRESSURE
cap_pres[ipore][j] = load_cap_pres(ipore, j, -1, sat_nodes_old[ipore][j]);
d_cap_pres_dS[ipore][j] = 0;
#else
cap_pres[ipore][j] = load_cap_pres(ipore, j, -1, sat_nodes[ipore][j]);
d_cap_pres_dS[ipore][j] = mp->d_cap_pres[var];
#endif
}
for (a = 0; a < DIM; a++) {
for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
95 changes: 89 additions & 6 deletions src/mm_std_models_shell.c
Original file line number Diff line number Diff line change
@@ -1959,6 +1959,7 @@ void porous_shell_open_source_model(
dbl mu = mp->viscosity; // Viscosity

dbl sat_nodes[MAX_POR_SHELL][MDE] = {{0.0}};

Check warning on line 1961 in src/mm_std_models_shell.c

GitHub Actions / build release unittest

variable ‘sat_nodes’ set but not used [-Wunused-but-set-variable]
dbl sat_nodes_old[MAX_POR_SHELL][MDE] = {{0.0}};
dbl cap_pres_nodes[MAX_POR_SHELL][MDE] = {{0.0}};
dbl d_cap_pres_nodes_dS[MAX_POR_SHELL][MDE] = {{0.0}};

@@ -1968,14 +1969,18 @@ void porous_shell_open_source_model(
dbl dHside_1_square_dS[MDE] = {0.0};
dbl sat_max_1 = 0.99;
dbl width_1 = 0.24;
// dbl width_1 = 0.0;
dbl alpha_1 = 0.5 * width_1;
dbl sat_center_1 = sat_max_1 - alpha_1;
dbl sat_normalized_1[MDE] = {0.0};

dbl Hside_2[MDE] = {0.0};
dbl dHside_2_dS[MDE] = {0.0};
dbl sat_min_2 = 0.5;
// dbl sat_min_2 = 0.1;
// dbl sat_min_2 = 5e-2;
dbl width_2 = 0.05;
// dbl width_2 = 0.00;
dbl alpha_2 = 0.5 * width_2;
dbl sat_center_2 = sat_min_2 - alpha_2;
dbl sat_normalized_2[MDE] = {0.0};
@@ -1986,10 +1991,14 @@ void porous_shell_open_source_model(
dbl dHside_3_square_dS[MDE] = {0.0};
dbl sat_max_3 = 0.99;
dbl width_3 = 0.24;
// dbl width_3 = 0.0;
dbl alpha_3 = 0.5 * width_3;
dbl sat_center_3 = sat_max_3 - alpha_3;
dbl sat_normalized_3[MDE] = {0.0};

#define LAG_SAT_NORM
#define LAG_CAP_PRESSURE

/* Extract all of the necessary information */
for (ipore = 0; ipore < pd->Num_Porous_Shell_Eqn; ipore++) {
var = porous_shell_var[ipore];
@@ -1998,19 +2007,39 @@ void porous_shell_open_source_model(
switch (ipore) {
case 0:
sat_nodes[ipore][j] = *esp->sh_sat_1[j];
sat_nodes_old[ipore][j] = *esp_old->sh_sat_1[j];
#ifdef LAG_SAT_NORM
sat_normalized_1[j] = sat_nodes_old[ipore][j] - sat_center_1;
#else
sat_normalized_1[j] = sat_nodes[ipore][j] - sat_center_1;
#endif
break;
case 1:
sat_nodes[ipore][j] = *esp->sh_sat_2[j];
sat_nodes_old[ipore][j] = *esp_old->sh_sat_2[j];
#ifdef LAG_SAT_NORM
sat_normalized_2[j] = sat_nodes_old[ipore][j] - sat_center_2;
#else
sat_normalized_2[j] = sat_nodes[ipore][j] - sat_center_2;
#endif
break;
case 2:
sat_nodes[ipore][j] = *esp->sh_sat_3[j];
sat_nodes_old[ipore][j] = *esp_old->sh_sat_3[j];
#ifdef LAG_SAT_NORM
sat_normalized_3[j] = sat_nodes_old[ipore][j] - sat_center_3;
#else
sat_normalized_3[j] = sat_nodes[ipore][j] - sat_center_3;
#endif
break;
}
#ifdef LAG_CAP_PRESSURE
cap_pres_nodes[ipore][j] = load_cap_pres(ipore, j, -1, sat_nodes_old[ipore][j]);
d_cap_pres_nodes_dS[ipore][j] = 0;
#else
cap_pres_nodes[ipore][j] = load_cap_pres(ipore, j, -1, sat_nodes[ipore][j]);
d_cap_pres_nodes_dS[ipore][j] = mp->d_cap_pres[var];
#endif
}
H[ipore] = porous_shell_height_model(ipore);
kappa[ipore] = porous_shell_cross_perm_model(ipore);
@@ -2020,12 +2049,22 @@ void porous_shell_open_source_model(
/* Apply heaviside function to deactivate flux at S_1 > 0.99*/

for (j = 0; j < ei[pg->imtrx]->dof[SHELL_SAT_1]; j++) {
if (sat_nodes[0][j] >= sat_max_1) {

// #define DISABLE_H1
// #define DISABLE_H2
// #define DISABLE_H3
#ifdef DISABLE_H1
Hside_1[j] = 1.0;
dHside_1_dS[j] = 0.0;
Hside_1_square[j] = 1.0;
dHside_1_square_dS[j] = 0.0;
#else
if (sat_nodes_old[0][j] >= sat_max_1) {
Hside_1[j] = 0.0;
dHside_1_dS[j] = 0.0;
Hside_1_square[j] = 0.0;
dHside_1_square_dS[j] = 0.0;
} else if (sat_nodes[0][j] <= (sat_max_1 - width_1)) {
} else if (sat_nodes_old[0][j] <= (sat_max_1 - width_1)) {
Hside_1[j] = 1.0;
dHside_1_dS[j] = 0.0;
Hside_1_square[j] = 1.0;
@@ -2035,35 +2074,53 @@ void porous_shell_open_source_model(
sin(M_PIE * sat_normalized_1[j] / alpha_1) / M_PIE);
dHside_1_dS[j] =
-0.5 * (1.0 / alpha_1 + cos(M_PIE * sat_normalized_1[j] / alpha_1) / alpha_1);
#ifdef LAG_SAT_NORM
dHside_1_dS[j] = 0;
#endif
Hside_1_square[j] = Hside_1[j] * Hside_1[j];
dHside_1_square_dS[j] = 2.0 * Hside_1[j] * dHside_1_dS[j];
}
#endif
}

/* Apply heaviside function to activate flux at S_2 > S_min*/
for (j = 0; j < ei[pg->imtrx]->dof[SHELL_SAT_2]; j++) {
if (sat_nodes[1][j] >= sat_min_2) {
#ifdef DISABLE_H2
Hside_2[j] = 1.0;
dHside_2_dS[j] = 0.0;
#else
if (sat_nodes_old[1][j] >= sat_min_2) {
Hside_2[j] = 1.0;
dHside_2_dS[j] = 0.0;
} else if (sat_nodes[1][j] <= (sat_min_2 - width_2)) {
} else if (sat_nodes_old[1][j] <= (sat_min_2 - width_2)) {
Hside_2[j] = 0.0;
dHside_2_dS[j] = 0.0;
} else {
Hside_2[j] = 0.5 * (1. + sat_normalized_2[j] / alpha_2 +
sin(M_PIE * sat_normalized_2[j] / alpha_2) / M_PIE);
dHside_2_dS[j] = 0.5 * (1.0 / alpha_2 + cos(M_PIE * sat_normalized_2[j] / alpha_2) / alpha_2);
#ifdef LAG_SAT_NORM
dHside_2_dS[j] = 0;
#endif
}
#endif
}

/* Apply heaviside function to deactivate flux at S_3 > 0.99*/
if (pd->e[pg->imtrx][R_SHELL_SAT_3]) {
for (j = 0; j < ei[pg->imtrx]->dof[SHELL_SAT_3]; j++) {
if (sat_nodes[2][j] >= sat_max_3) {
#ifdef DISABLE_H3
Hside_3[j] = 1.0;
dHside_3_dS[j] = 0.0;
Hside_3_square[j] = 1.0;
dHside_3_square_dS[j] = 0.0;
#else
if (sat_nodes_old[2][j] >= sat_max_3) {
Hside_3[j] = 0.0;
dHside_3_dS[j] = 0.0;
Hside_3_square[j] = 0.0;
dHside_3_square_dS[j] = 0.0;
} else if (sat_nodes[2][j] <= (sat_max_3 - width_3)) {
} else if (sat_nodes_old[2][j] <= (sat_max_3 - width_3)) {
Hside_3[j] = 1.0;
dHside_3_dS[j] = 0.0;
Hside_3_square[j] = 1.0;
@@ -2073,22 +2130,48 @@ void porous_shell_open_source_model(
sin(M_PIE * sat_normalized_3[j] / alpha_3) / M_PIE);
dHside_3_dS[j] =
-0.5 * (1.0 / alpha_3 + cos(M_PIE * sat_normalized_3[j] / alpha_3) / alpha_3);
#ifdef LAG_SAT_NORM
dHside_3_dS[j] = 0;
#endif
Hside_3_square[j] = Hside_3[j] * Hside_3[j];
dHside_3_square_dS[j] = 2.0 * Hside_3[j] * dHside_3_dS[j];
}
#endif
}
}

// #define DIODE_1
// #define DIODE_2
/* Populate the interporous flux */
for (j = 0; j < ei[pg->imtrx]->dof[SHELL_SAT_1]; j++) {
j_1_2[j] = (kappa[0] / mu) * (cap_pres_nodes[0][j] - cap_pres_nodes[1][j]) / (2.0 * H[0]);
j_1_2[j] += (kappa[1] / mu) * (cap_pres_nodes[0][j] - cap_pres_nodes[1][j]) / (2.0 * H[1]);
#ifdef DIODE_1
if (j_1_2[j] > 0.0) {
j_1_2[j] *= Hside_1_square[j] * Hside_2[j];
}
#elif defined DIODE_2
if (j_1_2[j] > 0.0) {
j_1_2[j] = 0;
}
#else
j_1_2[j] *= Hside_1_square[j] * Hside_2[j];
#endif

if (pd->e[pg->imtrx][R_SHELL_SAT_3]) {
j_2_3[j] = (kappa[1] / mu) * (cap_pres_nodes[2][j] - cap_pres_nodes[1][j]) / (2.0 * H[1]);
j_2_3[j] += (kappa[2] / mu) * (cap_pres_nodes[2][j] - cap_pres_nodes[1][j]) / (2.0 * H[2]);
#ifdef DIODE_1
if (j_2_3[j] > 0.0) {
j_2_3[j] *= Hside_2[j] * Hside_3_square[j];
}
#elif defined DIODE_2
if (j_2_3[j] > 0.0) {
j_2_3[j] = 0;
}
#else
j_2_3[j] *= Hside_2[j] * Hside_3_square[j];
#endif
}
}

0 comments on commit c438c17

Please sign in to comment.