From f2f7d6809442eddf1525ec00d2a73c3c1840d93f Mon Sep 17 00:00:00 2001 From: Ian Bertolacci Date: Tue, 18 Jun 2019 16:10:48 -0700 Subject: [PATCH] enable loop parallelism in NL_Function_Eval Loop calculating contributions from 2nd order derivatives and gravity had a scatter pattern that prevented parallelism. This is changed by: 1. storing the temporary values (u_right, u_front, u_upper) 2. Splitting scatter portion into a second loop below 3. reversing the accesses to perform a gather instead of a scatter --- pfsimulator/parflow_lib/nl_function_eval.c | 96 +++++++++++++++++----- 1 file changed, 77 insertions(+), 19 deletions(-) diff --git a/pfsimulator/parflow_lib/nl_function_eval.c b/pfsimulator/parflow_lib/nl_function_eval.c index ba64224ce..bd15b6212 100644 --- a/pfsimulator/parflow_lib/nl_function_eval.c +++ b/pfsimulator/parflow_lib/nl_function_eval.c @@ -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; @@ -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; @@ -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 */ @@ -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); @@ -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; @@ -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); @@ -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, @@ -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], @@ -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, @@ -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], @@ -721,7 +738,7 @@ 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], @@ -729,14 +746,54 @@ void NlFunctionEval(Vector * pressure, /* Current pressure values */ / 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] ); }); } @@ -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; @@ -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; } @@ -1732,6 +1793,3 @@ int NlFunctionEvalSizeOfTempData() { return 0; } - - -