[cig-commits] commit: Add a slope limiter for the pressure and viscosity*dv terms. It
Mercurial
hg at geodynamics.org
Sat Mar 12 05:36:59 PST 2011
changeset: 109:4c50d7d566ce
user: Walter Landry <wlandry at caltech.edu>
date: Tue Mar 08 20:14:39 2011 -0800
files: StokesFACOps.h StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C
description:
Add a slope limiter for the pressure and viscosity*dv terms. It
does not work, so it is turned off.
diff -r bb0d3477e3a5 -r 4c50d7d566ce StokesFACOps.h
--- a/StokesFACOps.h Sat Mar 05 23:47:30 2011 -0800
+++ b/StokesFACOps.h Tue Mar 08 20:14:39 2011 -0800
@@ -631,11 +631,60 @@ private:
return result;
}
+ /* A flux limiter for the pressure and viscosity derivatives. Minmod
+ for now.*/
+
+ double slope_limiter(const double &lower, const double &middle,
+ const double &upper)
+ {
+ if(lower>0 && middle>0 && upper>0)
+ return std::min(lower,std::min(middle,upper));
+ else if(lower<0 && middle<0 && upper<0)
+ return std::max(lower,std::max(middle,upper));
+ else return 0;
+ }
+
+ double dtau_xx_dx(pdat::SideData<double> &v,
+ pdat::CellData<double> &cell_viscosity,
+ const pdat::CellIndex ¢er,
+ const pdat::CellIndex &left,
+ const pdat::SideIndex ¢er_x,
+ const pdat::SideIndex &right_x,
+ const pdat::SideIndex &left_x,
+ const double &dx)
+ {
+ return 2*((v(right_x)-v(center_x))*cell_viscosity(center)
+ - (v(center_x)-v(left_x))*cell_viscosity(left))/(dx*dx);
+ }
+
+ double dtau_xy_dy(pdat::SideData<double> &v,
+ pdat::NodeData<double> &edge_viscosity,
+ const pdat::SideIndex ¢er_x,
+ const pdat::SideIndex &up_x,
+ const pdat::SideIndex &down_x,
+ const pdat::SideIndex ¢er_y,
+ const pdat::SideIndex ¢er_y_ip,
+ const pdat::SideIndex &up_y,
+ const pdat::SideIndex &up_y_ip,
+ const pdat::NodeIndex ¢er_e,
+ const pdat::NodeIndex &up_e,
+ const double &dx,
+ const double &dy)
+ {
+ return edge_viscosity(up_e)*((v(up_x)-v(center_x))/(dy*dy)
+ + (v(up_y)-v(up_y_ip))/(dx*dy))
+ - edge_viscosity(center_e)*((v(center_x)-v(down_x))/(dy*dy)
+ + (v(center_y)-v(center_y_ip))/(dx*dy));
+ }
+
/* The action of the velocity operator. It is written from the
perspective of vx, but pass in different values for center_x
etc. to get vy. */
- double v_operator(pdat::SideData<double> &v,
+ double v_operator(const hier::Box &pbox,
+ const int &axis,
+ const int &off_axis,
+ pdat::SideData<double> &v,
pdat::CellData<double> &p,
pdat::CellData<double> &cell_viscosity,
pdat::NodeData<double> &edge_viscosity,
@@ -651,21 +700,117 @@ private:
const pdat::NodeIndex ¢er_e,
const pdat::NodeIndex &up_e,
const hier::Index &ip,
+ const hier::Index &jp,
const double &dx,
const double &dy)
{
- const double dtau_xx_dx =
- 2*((v(right_x)-v(center_x))*cell_viscosity(center)
- - (v(center_x)-v(left_x))*cell_viscosity(left))/(dx*dx);
+ double dtau_xx_dx_lower, dtau_xx_dx_middle, dtau_xx_dx_upper,
+ dp_dx_lower, dp_dx_middle, dp_dx_upper;
+ dtau_xx_dx_middle=dtau_xx_dx(v,cell_viscosity,center,left,
+ center_x,right_x,left_x,dx);
+ dp_dx_middle=(p(center)-p(left))/dx;
+ if(center_x[axis]<=pbox.lower(axis)+1)
+ {
+ dtau_xx_dx_lower=dtau_xx_dx_middle;
+ dp_dx_lower=dp_dx_middle;
+ }
+ else
+ {
+ dtau_xx_dx_lower=dtau_xx_dx(v,cell_viscosity,center-ip,left-ip,
+ center_x-ip,right_x-ip,left_x-ip,dx);
+ dp_dx_lower=(p(center-ip)-p(left-ip))/dx;
+ }
+ if(center_x[axis]>=pbox.upper(axis))
+ {
+ dtau_xx_dx_upper=dtau_xx_dx_middle;
+ dp_dx_upper=dp_dx_middle;
+ }
+ else
+ {
+ dtau_xx_dx_upper=dtau_xx_dx(v,cell_viscosity,center+ip,left+ip,
+ center_x+ip,right_x+ip,left_x+ip,dx);
+ dp_dx_lower=(p(center+ip)-p(left+ip))/dx;
+ }
- const double dtau_xy_dy =
- edge_viscosity(up_e)*((v(up_x)-v(center_x))/(dy*dy)
- + (v(up_y)-v(up_y-ip))/(dx*dy))
- - edge_viscosity(center_e)*((v(center_x)-v(down_x))/(dy*dy)
- + (v(center_y)-v(center_y-ip))/(dx*dy));
- const double dp_dx=(p(center)-p(left))/dx;
+ const double dtau_xx_dx_limited(slope_limiter(dtau_xx_dx_lower,
+ dtau_xx_dx_middle,
+ dtau_xx_dx_upper));
+ const double
+ dp_dx_limited=slope_limiter(dp_dx_lower,dp_dx_middle,dp_dx_upper);
+
+
+ double dtau_xy_dy_lower, dtau_xy_dy_middle, dtau_xy_dy_upper;
+ dtau_xy_dy_middle=dtau_xy_dy(v,edge_viscosity,center_x,up_x,down_x,
+ center_y,center_y-ip,up_y,up_y-ip,
+ center_e,up_e,dx,dy);
+ if(center_y[off_axis]<=pbox.lower(off_axis)+1)
+ {
+ dtau_xy_dy_lower=dtau_xy_dy_middle;
+ }
+ else
+ {
+ dtau_xy_dy_lower=dtau_xy_dy(v,edge_viscosity,center_x-jp,up_x-jp,
+ down_x-jp,
+ center_y-jp,center_y-ip-jp,
+ up_y-jp,up_y-ip-jp,
+ center_e-jp,up_e-jp,dx,dy);
+ }
+ if(center_y[off_axis]>=pbox.lower(off_axis))
+ {
+ dtau_xy_dy_upper=dtau_xy_dy_middle;
+ }
+ else
+ {
+ dtau_xy_dy_upper=dtau_xy_dy(v,edge_viscosity,center_x+jp,up_x+jp,
+ down_x+jp,
+ center_y+jp,center_y-ip+jp,
+ up_y+jp,up_y-ip+jp,
+ up_e+jp,center_e+jp,dx,dy);
+ }
+
+ const double dtau_xy_dy_limited(slope_limiter(dtau_xy_dy_lower,
+ dtau_xy_dy_middle,
+ dtau_xy_dy_upper));
+
+ // return dtau_xx_dx_limited + dtau_xy_dy_limited - dp_dx_limited;
+
+ if(axis==0 && center[1]<=2)
+ {
+ tbox::plog << "v_op "
+ << center[0] << " "
+ << center[1] << " "
+ << axis << " "
+ << v(center_x) << " "
+ << v(center_y) << " "
+ << p(center) << " "
+ << (dtau_xx_dx_middle + dtau_xy_dy_middle - dp_dx_middle)
+ << " "
+ << (v(center_x + ip) - v(center_x)
+ + v(center_y+jp) - v(center_y)) << " "
+ // << dtau_xx_dx_middle << " "
+ // << dtau_xx_dx_limited << " "
+ // << dtau_xy_dy_middle << " "
+ // << dtau_xy_dy_limited << " "
+ // << dp_dx_middle << " "
+ // << dp_dx_limited << " "
+ << "\n";
+ }
+
+ return dtau_xx_dx_middle + dtau_xy_dy_middle - dp_dx_middle;
- return dtau_xx_dx + dtau_xy_dy - dp_dx;
+
+ // const double dtau_xx_dx =
+ // 2*((v(right_x)-v(center_x))*cell_viscosity(center)
+ // - (v(center_x)-v(left_x))*cell_viscosity(left))/(dx*dx);
+
+ // const double dtau_xy_dy =
+ // edge_viscosity(up_e)*((v(up_x)-v(center_x))/(dy*dy)
+ // + (v(up_y)-v(up_y-ip))/(dx*dy))
+ // - edge_viscosity(center_e)*((v(center_x)-v(down_x))/(dy*dy)
+ // + (v(center_y)-v(center_y-ip))/(dx*dy));
+ // const double dp_dx=(p(center)-p(left))/dx;
+
+ // return dtau_xx_dx + dtau_xy_dy - dp_dx;
}
/*!
diff -r bb0d3477e3a5 -r 4c50d7d566ce StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Sat Mar 05 23:47:30 2011 -0800
+++ b/StokesFACOps/Update_V.C Tue Mar 08 20:14:39 2011 -0800
@@ -68,8 +68,9 @@ void SAMRAI::solv::StokesFACOps::Update_
const double &theta_momentum)
{
const int off_axis=(axis==0) ? 1 : 0;
- hier::Index ip(0,0);
+ hier::Index ip(0,0), jp(0,0);
ip[axis]=1;
+ jp[off_axis]=1;
const pdat::SideIndex
center_x(center,axis,pdat::SideIndex::Lower),
@@ -117,9 +118,10 @@ void SAMRAI::solv::StokesFACOps::Update_
dx,dy);
double delta_Rx=v_rhs(center_x)
- - v_operator(v,p,cell_viscosity,edge_viscosity,center,left,center_x,
+ - v_operator(pbox,axis,off_axis,
+ v,p,cell_viscosity,edge_viscosity,center,left,center_x,
right_x,left_x,up_x,down_x,center_y,up_y,center_e,up_e,
- ip,dx,dy);
+ ip,jp,dx,dy);
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_Rx));
diff -r bb0d3477e3a5 -r 4c50d7d566ce StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Sat Mar 05 23:47:30 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Tue Mar 08 20:14:39 2011 -0800
@@ -254,9 +254,10 @@ void SAMRAI::solv::StokesFACOps::compute
else
{
v_resid(center_x)=v_rhs(center_x)
- - v_operator(v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator(pbox,0,1,
+ v,p,cell_viscosity,edge_viscosity,center,
left,center_x,right_x,left_x,up_x,down_x,
- center_y,up_y,center_e,up_e,ip,dx,dy);
+ center_y,up_y,center_e,up_e,ip,jp,dx,dy);
// tbox::plog << "vx "
// << v_resid(center_x) << " "
// << v(center_x) << " "
@@ -277,9 +278,11 @@ void SAMRAI::solv::StokesFACOps::compute
else
{
v_resid(center_y)=v_rhs(center_y)
- - v_operator(v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator(pbox,1,0,
+ v,p,cell_viscosity,edge_viscosity,center,
down,center_y,up_y,down_y,right_y,left_y,
- center_x,right_x,center_e,right_e,jp,dy,dx);
+ center_x,right_x,center_e,right_e,jp,ip,
+ dy,dx);
}
// tbox::plog << "vy "
// << v_resid(center_y) << " ";
More information about the CIG-COMMITS
mailing list