[cig-commits] commit: Remove the slope limiter, clean up things a bit and move dRc_dp and
Mercurial
hg at geodynamics.org
Sat Mar 12 05:37:05 PST 2011
changeset: 113:576b39c7e446
user: Walter Landry <wlandry at caltech.edu>
date: Thu Mar 10 20:44:11 2011 -0800
files: StokesFACOps.h StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/smooth_Tackley.C dRc_dp.h dRm_dv.h
description:
Remove the slope limiter, clean up things a bit and move dRc_dp and
dRm_dv into their own files.
diff -r 1d049d0f2567 -r 576b39c7e446 StokesFACOps.h
--- a/StokesFACOps.h Thu Mar 10 16:20:02 2011 -0800
+++ b/StokesFACOps.h Thu Mar 10 20:44:11 2011 -0800
@@ -563,128 +563,11 @@ private:
pdat::NodeData<double> &edge_viscosity,
const double &theta_momentum);
- /* The derivative of the momentum equation w/respect to velocity. It
- is written from the perspective of vx(center_x), but pass in
- different values for center etc. to get vy or vx(!center_x). */
-
- double dRm_dv(pdat::CellData<double> &cell_viscosity,
- pdat::NodeData<double> &edge_viscosity,
- const pdat::CellIndex ¢er,
- const pdat::CellIndex &left,
- const pdat::NodeIndex &up_e,
- const pdat::NodeIndex ¢er_e,
- const double &dx,
- const double &dy)
- {
- return -2*(cell_viscosity(center) + cell_viscosity(left))/(dx*dx)
- - (edge_viscosity(up_e) + edge_viscosity(center_e))/(dy*dy);
- }
-
- /* The derivative of the continuity equation with respect to
- pressure. Note that pressure does not appear in the continuity
- equation, so we use Tackley's method to chain together
- derivatives */
-
- 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;
- }
-
- /* 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(const hier::Box &pbox,
- const int &axis,
- const int &off_axis,
- pdat::SideData<double> &v,
+ double v_operator(pdat::SideData<double> &v,
pdat::CellData<double> &p,
pdat::CellData<double> &cell_viscosity,
pdat::NodeData<double> &edge_viscosity,
@@ -700,117 +583,19 @@ 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)
{
- 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;
- }
+ 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_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));
+ 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;
-
-
- // 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;
+ return dtau_xx_dx + dtau_xy_dy - dp_dx;
}
/*!
diff -r 1d049d0f2567 -r 576b39c7e446 StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Thu Mar 10 16:20:02 2011 -0800
+++ b/StokesFACOps/Update_V.C Thu Mar 10 20:44:11 2011 -0800
@@ -41,6 +41,7 @@
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
#include "Boundary.h"
+#include "dRc_dp.h"
/*
********************************************************************
* Updates one component of the velocity during a red-black *
@@ -118,10 +119,9 @@ void SAMRAI::solv::StokesFACOps::Update_
dx,dy);
double delta_Rx=v_rhs(center_x)
- - v_operator(pbox,axis,off_axis,
- v,p,cell_viscosity,edge_viscosity,center,left,center_x,
+ - v_operator(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,jp,dx,dy);
+ ip,dx,dy);
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_Rx));
diff -r 1d049d0f2567 -r 576b39c7e446 StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Thu Mar 10 16:20:02 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Thu Mar 10 20:44:11 2011 -0800
@@ -41,6 +41,7 @@
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
#include "Boundary.h"
+#include "dRc_dp.h"
/*
********************************************************************
* FACOperatorStrategy virtual *
@@ -240,10 +241,9 @@ void SAMRAI::solv::StokesFACOps::compute
else
{
v_resid(center_x)=v_rhs(center_x)
- - v_operator(pbox,0,1,
- v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator(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,jp,dx,dy);
+ center_y,up_y,center_e,up_e,ip,dx,dy);
// tbox::plog << "vx "
// << v_resid(center_x) << " "
// << v(center_x) << " "
@@ -264,10 +264,9 @@ void SAMRAI::solv::StokesFACOps::compute
else
{
v_resid(center_y)=v_rhs(center_y)
- - v_operator(pbox,1,0,
- v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator(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,ip,
+ center_x,right_x,center_e,right_e,jp,
dy,dx);
}
// tbox::plog << "vy "
diff -r 1d049d0f2567 -r 576b39c7e446 StokesFACOps/smooth_Tackley.C
--- a/StokesFACOps/smooth_Tackley.C Thu Mar 10 16:20:02 2011 -0800
+++ b/StokesFACOps/smooth_Tackley.C Thu Mar 10 20:44:11 2011 -0800
@@ -41,6 +41,7 @@
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
#include "Boundary.h"
+#include "dRc_dp.h"
/*
********************************************************************
* Workhorse function to smooth error using red-black *
diff -r 1d049d0f2567 -r 576b39c7e446 dRc_dp.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/dRc_dp.h Thu Mar 10 20:44:11 2011 -0800
@@ -0,0 +1,62 @@
+#ifndef GAMR_DRC_DP
+#define GAMR_DRC_DP
+
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/NodeData.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "dRm_dv.h"
+
+/* The derivative of the continuity equation with respect to
+ pressure. Note that pressure does not appear in the continuity
+ equation, so we use Tackley's method to chain together
+ derivatives */
+
+inline double dRc_dp(const SAMRAI::hier::Box &pbox,
+ const SAMRAI::pdat::CellIndex ¢er,
+ const SAMRAI::pdat::CellIndex &left,
+ const SAMRAI::pdat::CellIndex &right,
+ const SAMRAI::pdat::CellIndex &down,
+ const SAMRAI::pdat::CellIndex &up,
+ const SAMRAI::pdat::SideIndex &left_x,
+ const SAMRAI::pdat::SideIndex &right_x,
+ const SAMRAI::pdat::SideIndex &down_y,
+ const SAMRAI::pdat::SideIndex &up_y,
+ SAMRAI::pdat::CellData<double> &cell_viscosity,
+ SAMRAI::pdat::NodeData<double> &edge_viscosity,
+ SAMRAI::pdat::SideData<double> &v,
+ const double &dx,
+ const double &dy)
+{
+ const SAMRAI::pdat::NodeIndex
+ center_e(center,SAMRAI::pdat::NodeIndex::LowerLeft),
+ up_e(up,SAMRAI::pdat::NodeIndex::LowerLeft),
+ right_e(right,SAMRAI::pdat::NodeIndex::LowerLeft);
+ const SAMRAI::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;
+}
+
+
+#endif
diff -r 1d049d0f2567 -r 576b39c7e446 dRm_dv.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/dRm_dv.h Thu Mar 10 20:44:11 2011 -0800
@@ -0,0 +1,27 @@
+#ifndef GAMR_DRM_DV
+#define GAMR_DRM_DV
+
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/NodeData.h"
+
+/* The derivative of the momentum equation w/respect to velocity. It
+ is written from the perspective of vx(center_x), but pass in
+ different values for center etc. to get vy or vx(!center_x). */
+
+inline double dRm_dv(SAMRAI::pdat::CellData<double> &cell_viscosity,
+ SAMRAI::pdat::NodeData<double> &edge_viscosity,
+ const SAMRAI::pdat::CellIndex ¢er,
+ const SAMRAI::pdat::CellIndex &left,
+ const SAMRAI::pdat::NodeIndex &up_e,
+ const SAMRAI::pdat::NodeIndex ¢er_e,
+ const double &dx,
+ const double &dy)
+{
+ return -2*(cell_viscosity(center) + cell_viscosity(left))/(dx*dx)
+ - (edge_viscosity(up_e) + edge_viscosity(center_e))/(dy*dy);
+}
+
+
+
+
+#endif
More information about the CIG-COMMITS
mailing list