Skip to content
This repository was archived by the owner on May 20, 2024. It is now read-only.

enable loop parallelism in NL_Function_Eval #12

Closed
ian-bertolacci wants to merge 1 commit into from
Closed

enable loop parallelism in NL_Function_Eval #12

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
96 changes: 77 additions & 19 deletions pfsimulator/parflow_lib/nl_function_eval.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,11 +144,14 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */

/* Overland flow variables */ //sk
Vector *KW, *KE, *KN, *KS;
Vector *u_right_vec, *u_front_vec, *u_upper_vec; // @IJB
Vector *qx, *qy;
Subvector *kw_sub, *ke_sub, *kn_sub, *ks_sub, *qx_sub, *qy_sub;
Subvector *u_right_sub, *u_front_sub, *u_upper_sub;
Subvector *x_sl_sub, *y_sl_sub, *mann_sub;
Subvector *obf_sub;
double *kw_, *ke_, *kn_, *ks_, *qx_, *qy_;
double *u_right, *u_front, *u_upper; // @IJB
double *x_sl_dat, *y_sl_dat, *mann_dat;
double *obf_dat;
double q_overlnd;
Expand Down Expand Up @@ -198,12 +201,11 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
int nx_f, ny_f, nz_f;
int nx_p, ny_p, nz_p;
int nx_po, ny_po, nz_po;
int sy_p, sz_p;
int sx_p, sy_p, sz_p;
int ip, ipo, io;
int diffusive; //@RMM

double dtmp, dx, dy, dz, vol, ffx, ffy, ffz;
double u_right, u_front, u_upper;
double diff = 0.0e0;
double updir = 0.0e0;
double lower_cond, upper_cond;
Expand Down Expand Up @@ -246,6 +248,11 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
qx = NewVectorType(grid2d, 1, 1, vector_cell_centered_2D);
qy = NewVectorType(grid2d, 1, 1, vector_cell_centered_2D);

// @IJB
u_right_vec = NewVectorType(grid, 1, 1, fval->type);
u_front_vec = NewVectorType(grid, 1, 1, fval->type);
u_upper_vec = NewVectorType(grid, 1, 1, fval->type);


/* Calculate pressure dependent properties: density and saturation */

Expand Down Expand Up @@ -570,6 +577,11 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
x_ssl_sub = VectorSubvector(x_ssl, is);
y_ssl_sub = VectorSubvector(y_ssl, is);

/* @IJB added to enable parallelism */
u_right_sub = VectorSubvector(u_right_vec, is);
u_front_sub = VectorSubvector(u_front_vec, is);
u_upper_sub = VectorSubvector(u_upper_vec, is);

/* @RMM added to provide access to zmult */
z_mult_sub = VectorSubvector(z_mult, is);

Expand All @@ -596,6 +608,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
ny_p = SubvectorNY(p_sub);
nz_p = SubvectorNZ(p_sub);

sx_p = 1;
sy_p = nx_p;
sz_p = ny_p * nx_p;

Expand All @@ -616,12 +629,16 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
x_ssl_dat = SubvectorData(x_ssl_sub);
y_ssl_dat = SubvectorData(y_ssl_sub);

/* @IJB added to enable parallelism */
u_right = SubvectorData(u_right_sub);
u_front = SubvectorData(u_front_sub);
u_upper = SubvectorData(u_upper_sub);

/* @RMM added to provide variable dz */
z_mult_dat = SubvectorData(z_mult_sub);

qx_sub = VectorSubvector(qx, is);


GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx, ny, nz,
{
ip = SubvectorEltIndex(p_sub, i, j, k);
Expand Down Expand Up @@ -657,7 +674,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
diff = pp[ip] - pp[ip + 1];
updir = (diff / dx) * x_dir_g_c - x_dir_g;

u_right = z_mult_dat[ip] * ffx * del_y_slope * PMean(pp[ip], pp[ip + 1],
u_right[ip] = z_mult_dat[ip] * ffx * del_y_slope * PMean(pp[ip], pp[ip + 1],
permxp[ip], permxp[ip + 1])
* (diff / (dx * del_x_slope)) * x_dir_g_c
* RPMean(updir, 0.0,
Expand All @@ -670,7 +687,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
* upwind on pressure is currently implemented
* Sx < 0 implies flow goes left to right */

u_right += z_mult_dat[ip] * ffx * del_y_slope * PMean(pp[ip], pp[ip + 1],
u_right[ip] += z_mult_dat[ip] * ffx * del_y_slope * PMean(pp[ip], pp[ip + 1],
permxp[ip], permxp[ip + 1])
* (-x_dir_g)
* RPMean(updir, 0.0, rpp[ip] * dp[ip],
Expand All @@ -683,7 +700,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
diff = pp[ip] - pp[ip + sy_p];
updir = (diff / dy) * y_dir_g_c - y_dir_g;

u_front = z_mult_dat[ip] * ffy * del_x_slope
u_front[ip] = z_mult_dat[ip] * ffy * del_x_slope
* PMean(pp[ip], pp[ip + sy_p], permyp[ip], permyp[ip + sy_p])
* (diff / (dy * del_y_slope)) * y_dir_g_c
* RPMean(updir, 0.0,
Expand All @@ -697,7 +714,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
* Sy < 0 implies flow goes from left to right
*/

u_front += z_mult_dat[ip] * ffy * del_x_slope
u_front[ip] += z_mult_dat[ip] * ffy * del_x_slope
* PMean(pp[ip], pp[ip + sy_p], permyp[ip], permyp[ip + sy_p])
* (-y_dir_g)
* RPMean(updir, 0.0, rpp[ip] * dp[ip],
Expand All @@ -721,22 +738,62 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */

diff = (lower_cond - upper_cond);

u_upper = ffz * del_x_slope * del_y_slope
u_upper[ip] = ffz * del_x_slope * del_y_slope
* PMeanDZ(permzp[ip], permzp[ip + sz_p], z_mult_dat[ip], z_mult_dat[ip + sz_p])
* diff
* RPMean(lower_cond, upper_cond, rpp[ip] * dp[ip],
rpp[ip + sz_p] * dp[ip + sz_p])
/ viscosity;

/* velocity data jjb */
vx[vxi] = u_right / ffx;
vy[vyi] = u_front / ffy;
vz[vzi] = u_upper / ffz;

fp[ip] += dt * (u_right + u_front + u_upper);
fp[ip + 1] -= dt * u_right;
fp[ip + sy_p] -= dt * u_front;
fp[ip + sz_p] -= dt * u_upper;
vx[vxi] = u_right[ip] / ffx;
vy[vyi] = u_front[ip] / ffy;
vz[vzi] = u_upper[ip] / ffz;

fp[ip] += dt * (u_right[ip] + u_front[ip] + u_upper[ip]);
});
}

// Gather portion
ForSubgridI(is, GridSubgrids(grid))
{
subgrid = GridSubgrid(grid, is);

p_sub = VectorSubvector(pressure, is);
f_sub = VectorSubvector(fval, is);

/* @IJB added to enable parallelism */
u_right_sub = VectorSubvector(u_right_vec, is);
u_front_sub = VectorSubvector(u_front_vec, is);
u_upper_sub = VectorSubvector(u_upper_vec, is);

ix = SubgridIX(subgrid) - 1;
iy = SubgridIY(subgrid) - 1;
iz = SubgridIZ(subgrid) - 1;

nx = SubgridNX(subgrid) + 1;
ny = SubgridNY(subgrid) + 1;
nz = SubgridNZ(subgrid) + 1;

nx_p = SubvectorNX(p_sub);
ny_p = SubvectorNY(p_sub);
nz_p = SubvectorNZ(p_sub);

sx_p = 1;
sy_p = nx_p;
sz_p = ny_p * nx_p;

fp = SubvectorData(f_sub);

/* @IJB added to enable parallelism */
u_right = SubvectorData(u_right_sub);
u_front = SubvectorData(u_front_sub);
u_upper = SubvectorData(u_upper_sub);

GrGeomInLoop(i, j, k, gr_domain, r, ix, iy, iz, nx+1, ny+1, nz+1,
{
ip = SubvectorEltIndex(p_sub, i, j, k);
fp[ip] += dt*( -u_right[ip - sx_p] + -u_front[ip - sy_p] + -u_upper[ip - sz_p] );
});
}

Expand Down Expand Up @@ -802,6 +859,7 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
ny_p = SubvectorNY(p_sub);
nz_p = SubvectorNZ(p_sub);

sx_p = 1;
sy_p = nx_p;
sz_p = ny_p * nx_p;

Expand Down Expand Up @@ -1579,6 +1637,9 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */
FreeVector(KS);
FreeVector(qx);
FreeVector(qy);
FreeVector(u_right_vec);
FreeVector(u_front_vec);
FreeVector(u_upper_vec);

return;
}
Expand Down Expand Up @@ -1732,6 +1793,3 @@ int NlFunctionEvalSizeOfTempData()
{
return 0;
}