[cig-commits] commit: Fix a bug where dRm/dp was using pressures outside the domain. Also factor out dRm/dp into a separate function.
Mercurial
hg at geodynamics.org
Sat Feb 26 00:01:57 PST 2011
changeset: 101:bd8a792b4a43
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Feb 25 23:59:46 2011 -0800
files: StokesFACOps.h StokesFACOps/smooth_Tackley.C
description:
Fix a bug where dRm/dp was using pressures outside the domain. Also factor out dRm/dp into a separate function.
diff -r a97ac60f56f4 -r bd8a792b4a43 StokesFACOps.h
--- a/StokesFACOps.h Fri Feb 25 23:57:22 2011 -0800
+++ b/StokesFACOps.h Fri Feb 25 23:59:46 2011 -0800
@@ -581,6 +581,52 @@ private:
- (edge_viscosity(up_e) + edge_viscosity(center_e))/(dy*dy);
}
+ double dRc_dp(const hier::Box &pbox,
+ const pdat::CellIndex ¢er,
+ const pdat::CellIndex &left,
+ const pdat::CellIndex &right,
+ const pdat::CellIndex &down,
+ const pdat::CellIndex &up,
+ const pdat::SideIndex &left_x,
+ const pdat::SideIndex &right_x,
+ const pdat::SideIndex &down_y,
+ const pdat::SideIndex &up_y,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::NodeData<double> &edge_viscosity,
+ pdat::SideData<double> &v,
+ const double &dx,
+ const double &dy)
+ {
+ const pdat::NodeIndex center_e(center,pdat::NodeIndex::LowerLeft),
+ up_e(up,pdat::NodeIndex::LowerLeft),
+ right_e(right,pdat::NodeIndex::LowerLeft);
+ const hier::Index ip(1,0), jp(0,1);
+ const double dRm_dp_xp(1/dx), dRm_dp_xm(-1/dx),
+ dRm_dp_yp(1/dy), dRm_dp_ym(-1/dy),
+ dRc_dvx_p(-1/dx), dRc_dvx_m(1/dx),
+ dRc_dvy_p(-1/dy), dRc_dvy_m(1/dy);
+
+ double result(0);
+
+ if(!(center[0]==pbox.lower(0) && v(left_x)==boundary_value))
+ result+=dRc_dvx_p * dRm_dp_xp/dRm_dv(cell_viscosity,edge_viscosity,right,
+ center,up_e+ip,center_e+ip,dx,dy);
+
+ if(!(center[0]==pbox.upper(0)+1 && v(right_x)==boundary_value))
+ result+=dRc_dvx_m * dRm_dp_xm/dRm_dv(cell_viscosity,edge_viscosity,
+ center,left,up_e,center_e,dx,dy);
+
+ if(!(center[1]==pbox.lower(1) && v(down_y)==boundary_value))
+ result+=dRc_dvy_p * dRm_dp_yp/dRm_dv(cell_viscosity,edge_viscosity,up,
+ center,right_e+jp,center_e+jp,dy,dx);
+
+ if(!(center[1]==pbox.upper(1)+1 && v(up_y)==boundary_value))
+ result+=dRc_dvy_m * dRm_dp_ym/dRm_dv(cell_viscosity,edge_viscosity,
+ center,down,right_e,center_e,dy,dx);
+
+ return result;
+ }
+
/*!
* @brief Solve the coarsest level using HYPRE
*/
diff -r a97ac60f56f4 -r bd8a792b4a43 StokesFACOps/smooth_Tackley.C
--- a/StokesFACOps/smooth_Tackley.C Fri Feb 25 23:57:22 2011 -0800
+++ b/StokesFACOps/smooth_Tackley.C Fri Feb 25 23:59:46 2011 -0800
@@ -92,8 +92,6 @@ void SAMRAI::solv::StokesFACOps::smooth_
double theta_momentum=0.7;
double theta_continuity=1.0;
-
- const hier::Index ip(1,0), jp(0,1);
/*
* Smooth the number of sweeps specified or until
@@ -454,25 +452,20 @@ void SAMRAI::solv::StokesFACOps::smooth_
++right[0];
--left[0];
- const pdat::NodeIndex
- center_e(center,pdat::NodeIndex::LowerLeft),
- up_e(up,pdat::NodeIndex::LowerLeft),
- right_e(right,pdat::NodeIndex::LowerLeft);
+ const pdat::SideIndex
+ center_x(center,0,pdat::SideIndex::Lower),
+ left_x(left,0,pdat::SideIndex::Lower),
+ right_x(right,0,pdat::SideIndex::Lower),
+ center_y(center,1,pdat::SideIndex::Lower),
+ up_y(up,1,pdat::SideIndex::Lower),
+ down_y(down,1,pdat::SideIndex::Lower);
/* Update p */
if(set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))])
{
- double dvx_dx=
- (v(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Upper))
- - v(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))/dx;
- double dvy_dy=
- (v(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Upper))
- - v(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))/dy;
+ double dvx_dx=(v(right_x) - v(center_x))/dx;
+ double dvy_dy=(v(up_y) - v(center_y))/dy;
double delta_R_continuity=
p_rhs(center) - dvx_dx - dvy_dy;
@@ -480,39 +473,10 @@ void SAMRAI::solv::StokesFACOps::smooth_
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_R_continuity));
- const double dRm_dp_xp(1/dx), dRm_dp_xm(-1/dx),
- dRm_dp_yp(1/dy), dRm_dp_ym(-1/dy),
- dRc_dvx_p(-1/dx), dRc_dvx_m(1/dx),
- dRc_dvy_p(-1/dy), dRc_dvy_m(1/dy);
-
- const double dRm_dvx_p =
- dRm_dv(cell_viscosity,edge_viscosity,
- right,center,up_e+ip,center_e+ip,dx,dy);
-
- const double dRm_dvx_m =
- dRm_dv(cell_viscosity,edge_viscosity,
- center,left,up_e,center_e,dx,dy);
-
- const double dRm_dvy_p =
- dRm_dv(cell_viscosity,edge_viscosity,
- up,center,right_e+jp,center_e+jp,dy,dx);
-
- const double dRm_dvy_m =
- dRm_dv(cell_viscosity,edge_viscosity,
- center,down,right_e,center_e,dy,dx);
-
- const double dRc_dp=dRc_dvx_p * dRm_dp_xp/dRm_dvx_p
- + dRc_dvx_m * dRm_dp_xm/dRm_dvx_m
- + dRc_dvy_p * dRm_dp_yp/dRm_dvy_p
- + dRc_dvy_m * dRm_dp_ym/dRm_dvy_m;
-
-
- dp(center)=
- delta_R_continuity*theta_continuity/dRc_dp;
- // dp(center)=
- // delta_R_continuity*theta_continuity;
-
-
+ dp(center)=delta_R_continuity*theta_continuity
+ /dRc_dp(pbox,center,left,right,down,up,
+ left_x,right_x,down_y,up_y,
+ cell_viscosity,edge_viscosity,v,dx,dy);
// if(ln==2)
// tbox::plog << "smooth p "
More information about the CIG-COMMITS
mailing list