Skip to content

Commit

Permalink
standalone qtensor solve
Browse files Browse the repository at this point in the history
  • Loading branch information
@wortiz
wortiz committed Jan 16, 2024
1 parent 6b9706b commit 43f3b60
Show file tree
Hide file tree
Showing 17 changed files with 785 additions and 13 deletions.
9 changes: 9 additions & 0 deletions include/mm_as_structs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -488,6 +488,7 @@ struct Element_Variable_Pointers {
dbl *P_star[MDE]; /* pressure */
dbl *S[MAX_MODES][DIM][DIM][MDE]; /* polymeric stress tensor, for each mode */
dbl *G[DIM][DIM][MDE]; /* velocity gradient tensor */
dbl *Q[DIM][DIM][MDE]; /* Q Tensor */
dbl *F[MDE]; /* Fill */
dbl *V[MDE]; /* Potential; added by KSC: 2/3/99 */
dbl *qs[MDE]; /* Surface charge density */
Expand Down Expand Up @@ -610,6 +611,7 @@ struct Element_Stiffness_Pointers {
stress tensor */
dbl ****G; /* *G[DIM][DIM][MDE], velocity gradient
tensor */
dbl ****Q; /* *G[DIM][DIM][MDE], Q Tensor */
dbl **V; /* *V[MDE], voltage potential */
dbl **qs; /* *qs[MDE], surface charge density */
dbl **F; /* *F[MDE], fill */
Expand Down Expand Up @@ -1647,6 +1649,7 @@ struct Field_Variables {
dbl P_star;
dbl S[MAX_MODES][DIM][DIM]; /* Polymer Stress, for each mode */
dbl G[DIM][DIM]; /* Velocity Gradient */
dbl Q[DIM][DIM]; /* Q Tensor */
dbl F; /* Fill */
dbl V; /* Voltage */
dbl qs; /* Surface charge density (shell element) */
Expand Down Expand Up @@ -1832,6 +1835,8 @@ struct Field_Variables {
dbl grad_S[MAX_MODES][DIM][DIM][DIM]; /* Gradient of polymer stress tensor( or most of it!) */
dbl div_S[MAX_MODES][DIM]; /* Divergence of polymer stress tensor */
dbl grad_G[DIM][DIM][DIM]; /* Gradient of velocity tensor ( or most of it!) */
dbl grad_Q[DIM][DIM][DIM]; /* Gradient of Q tensor ( or most of it!) */
dbl div_Q[DIM]; /* divergence of Q tensor */
dbl grad_Gt[DIM][DIM][DIM]; /* Gradient of the transpose of the velocity tensor */
dbl div_G[DIM]; /* Divergence of velocity gradient tensor */
dbl div_Gt[DIM]; /* Divergence of the transpose of velocity gradient tensor */
Expand Down Expand Up @@ -1911,6 +1916,9 @@ struct Field_Variables {
dbl d_grad_S_dmesh[MAX_MODES][DIM][DIM][DIM][DIM][MDE];
dbl d_div_S_dmesh[MAX_MODES][DIM][DIM][MDE];

dbl d_grad_Q_dmesh[DIM][DIM][DIM][DIM][MDE];
dbl d_div_Q_dmesh[DIM][DIM][MDE];

dbl d_grad_G_dmesh[DIM][DIM][DIM][DIM][MDE];
dbl d_div_G_dmesh[DIM][DIM][MDE];
dbl d_div_Gt_dmesh[DIM][DIM][MDE];
Expand Down Expand Up @@ -2051,6 +2059,7 @@ struct Diet_Field_Variables {
dbl n[DIM]; /* normal vector to level set field OR shell normal */
dbl S[MAX_MODES][DIM][DIM]; /* Polymer Stress, for each modes */
dbl G[DIM][DIM]; /* Velocity Gradient */
dbl Q[DIM][DIM]; /* Q tensor */
dbl nn; /* This is the bond evolution */
dbl p_liq; /* porous media liq-pressure variable. */
dbl p_gas; /* porous media gas-pressure variable. */
Expand Down
24 changes: 24 additions & 0 deletions include/mm_names.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7501,6 +7501,12 @@ struct Equation_Names EQ_Name[] = {
{"R_VSTAR", "VSTAR", R_VSTAR},
{"R_WSTAR", "WSTAR", R_WSTAR},
{"R_EDDY_NU", "EDDY_NU", R_EDDY_NU},
{"R_QTENSOR11", "QTENSOR11", QTENSOR11},
{"R_QTENSOR12", "QTENSOR12", QTENSOR12},
{"R_QTENSOR13", "QTENSOR13", QTENSOR13},
{"R_QTENSOR22", "QTENSOR22", QTENSOR22},
{"R_QTENSOR23", "QTENSOR23", QTENSOR23},
{"R_QTENSOR33", "QTENSOR33", QTENSOR33},

/*
* Note -> these entries must remain until we get rid
Expand Down Expand Up @@ -7801,6 +7807,12 @@ struct Equation_Names Var_Name[] = {
{"VSTAR", "USY", VSTAR},
{"WSTAR", "USZ", WSTAR},
{"EDDY_NU", "EDDY_NU", EDDY_NU}, // 214
{"QTENSOR11", "Q11", QTENSOR11}, // 215
{"QTENSOR12", "Q12", QTENSOR12},
{"QTENSOR13", "Q13", QTENSOR13},
{"QTENSOR22", "Q22", QTENSOR22},
{"QTENSOR23", "Q23", QTENSOR23},
{"QTENSOR33", "Q33", QTENSOR33},

{"MESH_POSITION1", "X", MESH_POSITION1},
{"MESH_POSITION2", "Y", MESH_POSITION2}, /* 216 */
Expand Down Expand Up @@ -8072,6 +8084,12 @@ struct Equation_Names Exo_Var_Names[] = {
{"V Int.", "USY", VSTAR},
{"W Int.", "USZ", WSTAR},
{"Eddy Turbulence Viscosity.", "EDDY_NU", EDDY_NU},
{"Q Tensor 11", "QXX", QTENSOR11},
{"Q Tensor 12", "QXY", QTENSOR12},
{"Q Tensor 13", "QXZ", QTENSOR13},
{"Q Tensor 22", "QYY", QTENSOR22},
{"Q Tensor 23", "QYZ", QTENSOR23},
{"Q Tensor 33", "QZZ", QTENSOR33},
};

int Num_Exo_Var_Names = sizeof(Exo_Var_Names) / sizeof(struct Equation_Names);
Expand Down Expand Up @@ -8379,6 +8397,12 @@ struct Equation_Names Var_Units[] = {
{"VSTAR", "[1]", VSTAR},
{"WSTAR", "[1]", WSTAR},
{"EDDY_NU", "[1]", EDDY_NU},
{"QTENSOR11", "[1]", QTENSOR11},
{"QTENSOR12", "[1]", QTENSOR12},
{"QTENSOR13", "[1]", QTENSOR13},
{"QTENSOR22", "[1]", QTENSOR22},
{"QTENSOR23", "[1]", QTENSOR23},
{"QTENSOR33", "[1]", QTENSOR33},
};

int Num_Var_Units = sizeof(Var_Units) / sizeof(struct Equation_Names);
Expand Down
2 changes: 2 additions & 0 deletions include/mm_qtensor_model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,8 @@ EXTERN int bias_eigenvector_to(dbl *, /* eigenvector */

EXTERN void assemble_qtensor(dbl *);

EXTERN int assemble_qtensor_full_fill(void);

EXTERN void assemble_new_qtensor /* Ryan's qtensor */
(dbl *);

Expand Down
15 changes: 14 additions & 1 deletion include/rf_fem_const.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -519,6 +519,13 @@
#define VSTAR 212
#define WSTAR 213
#define EDDY_NU 214
#define QTENSOR11 215
#define QTENSOR12 216
#define QTENSOR13 217
#define QTENSOR22 218
#define QTENSOR23 219
#define QTENSOR33 220

/*
* define a variable to hold an external field which will be
* held fixed in the problem but parametered by the basis functions
Expand Down Expand Up @@ -925,7 +932,13 @@
#define R_VSTAR 212
#define R_WSTAR 213
#define R_EDDY_NU 214
#define V_LAST 215
#define R_QTENSOR11 215
#define R_QTENSOR12 216
#define R_QTENSOR13 217
#define R_QTENSOR22 218
#define R_QTENSOR23 219
#define R_QTENSOR33 220
#define V_LAST 221

/* MMH
* This is used for those parts of the code that want to ensure
Expand Down
30 changes: 30 additions & 0 deletions src/bc_colloc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2587,6 +2587,36 @@ int load_variable(double *x_var, /* variable value */
var = EDDY_NU;
*d_x_var = 1.;
break;
case QTENSOR11:
*x_var = fv->Q[0][0];
var = QTENSOR11;
*d_x_var = 1.;
break;
case QTENSOR12:
*x_var = fv->Q[0][1];
var = QTENSOR12;
*d_x_var = 1.;
break;
case QTENSOR13:
*x_var = fv->Q[0][2];
var = QTENSOR13;
*d_x_var = 1.;
break;
case QTENSOR22:
*x_var = fv->Q[1][1];
var = QTENSOR22;
*d_x_var = 1.;
break;
case QTENSOR23:
*x_var = fv->Q[1][2];
var = QTENSOR23;
*d_x_var = 1.;
break;
case QTENSOR33:
*x_var = fv->Q[2][2];
var = QTENSOR33;
*d_x_var = 1.;
break;
case LIGHT_INTP:
*x_var = fv->poynt[0];
var = LIGHT_INTP;
Expand Down
172 changes: 171 additions & 1 deletion src/load_field_variables.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include "mm_mp.h"
#include "mm_post_def.h"
#include "rf_fem.h"
#include "rf_fem_const.h"
#include "std.h"

/***************************************************************************/
Expand Down Expand Up @@ -1278,6 +1279,38 @@ int load_fv(void)
}
}

int R_Q[DIM][DIM];
R_Q[0][0] = QTENSOR11;
R_Q[0][1] = QTENSOR12;
R_Q[0][2] = QTENSOR13;
R_Q[1][0] = QTENSOR12;
R_Q[1][1] = QTENSOR22;
R_Q[1][2] = QTENSOR23;
R_Q[2][0] = QTENSOR13;
R_Q[2][1] = QTENSOR23;
R_Q[2][2] = QTENSOR33;
for (p = 0; pdgv[QTENSOR11] && p < VIM; p++) {
for (q = p; q < VIM; q++) {
v = R_Q[p][q];
if (pdgv[v]) {
fv->Q[p][q] = fv_old->Q[p][q] = fv_dot->Q[p][q] = 0.0;
if (p != q)
fv->Q[q][p] = fv_old->Q[q][p] = fv_dot->Q[q][p] = 0.0;
dofs = ei[upd->matrix_index[v]]->dof[v];
for (i = 0; i < dofs; i++) {
fv->Q[p][q] += *esp->Q[p][q][i] * bf[v]->phi[i];
if (pd->TimeIntegration != STEADY) {
fv_old->Q[p][q] += *esp_old->Q[p][q][i] * bf[v]->phi[i];
fv_dot->Q[p][q] += *esp_dot->Q[p][q][i] * bf[v]->phi[i];
}
}
fv->Q[q][p] = fv->Q[p][q];
fv_old->Q[q][p] = fv_old->Q[p][q];
fv_dot->Q[q][p] = fv_dot->Q[p][q];
}
}
}

/*
* Species Unknown Variable
* Modifications for non-dilute systems:
Expand Down Expand Up @@ -2605,7 +2638,7 @@ int load_fv_grads(void)
/*
* curl(v)
*/
if (CURL_V != -1 && !InShellElementWithParentElementCoverage) {
if (pd->gv[QTENSOR11] || (CURL_V != -1 && !InShellElementWithParentElementCoverage)) {
v = VELOCITY1;
bfn = bf[v];

Expand Down Expand Up @@ -2993,6 +3026,60 @@ int load_fv_grads(void)
}
}

if (pd->gv[QTENSOR11]) {
v = QTENSOR11;
dofs = ei[upd->matrix_index[v]]->dof[v];
for (p = 0; p < VIM; p++) {
for (q = 0; q < VIM; q++) {
for (r = 0; r < VIM; r++) {
fv->grad_Q[r][p][q] = 0.0;

for (i = 0; i < dofs; i++) {
fv->grad_Q[r][p][q] += *esp->Q[p][q][i] * bf[v]->grad_phi[i][r];
}
}
}
}

/*
* div(Q) - this is a vector!
*/
for (r = 0; r < dim; r++) {
fv->div_Q[r] = 0.0;
for (q = 0; q < dim; q++) {
fv->div_Q[r] += fv->grad_Q[q][q][r];
}
}

if (pd->CoordinateSystem != CARTESIAN) {
for (s = 0; s < VIM; s++) {
for (r = 0; r < VIM; r++) {
for (p = 0; p < VIM; p++) {
fv->div_Q[s] += fv->Q[p][s] * fv->grad_e[p][r][s];
}
}
}

for (s = 0; s < VIM; s++) {
for (r = 0; r < VIM; r++) {
for (q = 0; q < VIM; q++) {
fv->div_Q[s] += fv->Q[r][q] * fv->grad_e[q][r][s];
}
}
}
}

} else if (zero_unused_grads && upd->vp[pg->imtrx][QTENSOR11] == -1) {
for (p = 0; p < VIM; p++) {
fv->div_Q[p] = 0.0;
for (q = 0; q < VIM; q++) {
for (r = 0; r < VIM; r++) {
fv->grad_Q[r][p][q] = 0.;
}
}
}
}

/*
* grad(vd)
*/
Expand Down Expand Up @@ -6778,6 +6865,89 @@ int load_fv_mesh_derivs(int okToZero)
memset(fv->d_div_G_dmesh, 0, siz);
memset(fv->d_div_Gt_dmesh, 0, siz);
}
/*
* d(grad(Q))/dmesh
*/

if (pd->gv[QTENSOR11]) {
v = QTENSOR11;
bfv = bf[v];
vdofs = ei[upd->matrix_index[v]]->dof[v];

siz = sizeof(double) * DIM * DIM * DIM * DIM * MDE;
memset(fv->d_grad_Q_dmesh, 0, siz);
siz = sizeof(double) * DIM * DIM * MDE;
memset(fv->d_div_Q_dmesh, 0, siz);

for (p = 0; p < VIM; p++) {
for (q = 0; q < VIM; q++) {
for (r = 0; r < VIM; r++) {
for (b = 0; b < VIM; b++) {
for (j = 0; j < mdofs; j++) {
for (i = 0; i < vdofs; i++) {
fv->d_grad_Q_dmesh[r][p][q][b][j] +=
*esp->Q[p][q][i] * bfv->d_grad_phi_dmesh[i][r][b][j];
}
}
}
}
}
}

for (r = 0; r < VIM; r++) {
for (q = 0; q < VIM; q++) {
for (b = 0; b < VIM; b++) {
for (j = 0; j < mdofs; j++) {
fv->d_div_Q_dmesh[r][b][j] += fv->d_grad_Q_dmesh[q][q][r][b][j];
}
}
}
}

if (pd->CoordinateSystem != CARTESIAN) {
for (s = 0; s < VIM; s++) {
for (r = 0; r < VIM; r++) {
for (p = 0; p < VIM; p++) {
for (q = 0; q < VIM; q++) {
for (b = 0; b < dim; b++) {
for (j = 0; j < mdofs; j++) {
fv->d_div_Q_dmesh[s][b][j] +=
fv->Q[p][q] *
(fv->d_grad_e_dq[p][r][q][b] * bfm->phi[j] * (double)delta(s, q) +
fv->d_grad_e_dq[q][p][s][b] * bfm->phi[j] * (double)delta(r, p));
}
}
}
}
}
}
}

if (pd->CoordinateSystem != CARTESIAN) {
for (s = 0; s < VIM; s++) {
for (r = 0; r < VIM; r++) {
for (p = 0; p < VIM; p++) {
for (q = 0; q < VIM; q++) {
for (b = 0; b < VIM; b++) {
for (j = 0; j < mdofs; j++) {
fv->d_div_Q_dmesh[s][b][j] +=
fv->Q[q][p] *
(fv->d_grad_e_dq[p][r][q][b] * bfm->phi[j] * (double)delta(s, q) +
fv->d_grad_e_dq[q][p][s][b] * bfm->phi[j] * (double)delta(r, p));
}
}
}
}
}
}
}

} else if (upd->vp[pg->imtrx][QTENSOR11] != -1 && okToZero) {
siz = sizeof(double) * DIM * DIM * DIM * DIM * MDE;
memset(fv->d_grad_Q_dmesh, 0, siz);
siz = sizeof(double) * DIM * DIM * MDE;
memset(fv->d_div_Q_dmesh, 0, siz);
}

/*
* d(grad(vd))/dmesh
Expand Down
5 changes: 5 additions & 0 deletions src/mm_as_alloc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -943,6 +943,11 @@ int assembly_alloc(Exo_DB *exo)
(void)memset(esp->S[0][0][0], 0, MAX_MODES * VIM * VIM * MDE * sizeof(dbl *));
}

if (Num_Var_In_Type[imtrx][QTENSOR11]) {
esp->Q = (dbl ****)array_alloc(3, VIM, VIM, MDE, sizeof(dbl *));
(void)memset(esp->Q[0][0], 0, VIM * VIM * MDE * sizeof(dbl *));
}

/* VELOCITY_GRADIENT */
if (Num_Var_In_Type[imtrx][VELOCITY_GRADIENT11]) {
esp->G = (dbl ****)array_alloc(3, VIM, VIM, MDE, sizeof(dbl *));
Expand Down
Loading

0 comments on commit 43f3b60

Please sign in to comment.