[cig-commits] commit: Centralize the calculation of the velocity operator for the residual
Mercurial
hg at geodynamics.org
Sat Mar 5 20:38:59 PST 2011
changeset: 103:c82543f37f62
user: Walter Landry <wlandry at caltech.edu>
date: Sat Mar 05 11:29:15 2011 -0800
files: FACStokes/initializeLevelData.C StokesFACOps.h StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C
description:
Centralize the calculation of the velocity operator for the residual
and smoothing. This also fixes the residual for isoviscous and
variable viscosity. Also added code for a circular inclusion with a
different viscosity. For straight smoothing, no multigrid, this seems
to work well for viscosity ranges up to 10^10.
diff -r 087c7e707b4d -r c82543f37f62 FACStokes/initializeLevelData.C
--- a/FACStokes/initializeLevelData.C Sat Mar 05 11:25:10 2011 -0800
+++ b/FACStokes/initializeLevelData.C Sat Mar 05 11:29:15 2011 -0800
@@ -21,20 +21,7 @@
#include "SAMRAI/hier/Variable.h"
#include "SAMRAI/hier/VariableDatabase.h"
-extern "C" {
- void F77_FUNC(setexactandrhs2d, SETEXACTANDRHS2D) (const int & ifirst0,
- const int & ilast0,
- const int & ifirst1,
- const int & ilast1,
- double* exact,
- double* rhs,
- const double* dx,
- const double* xlower);
-}
-
-namespace SAMRAI {
-
- /*
+/*
*************************************************************************
* Initialize data on a level. *
* *
@@ -42,92 +29,164 @@ namespace SAMRAI {
* Fill the rhs and exact solution. *
*************************************************************************
*/
- void FACStokes::initializeLevelData
- (const tbox::Pointer<hier::BasePatchHierarchy> patch_hierarchy,
- const int level_number,
- const double init_data_time,
- const bool can_be_refined,
- const bool initial_time,
- const tbox::Pointer<hier::BasePatchLevel> old_level,
- const bool allocate_data)
- {
+void SAMRAI::FACStokes::initializeLevelData
+(const tbox::Pointer<hier::BasePatchHierarchy> patch_hierarchy,
+ const int level_number,
+ const double init_data_time,
+ const bool can_be_refined,
+ const bool initial_time,
+ const tbox::Pointer<hier::BasePatchLevel> old_level,
+ const bool allocate_data)
+{
- (void)init_data_time;
- (void)can_be_refined;
- (void)initial_time;
- (void)old_level;
+ (void)init_data_time;
+ (void)can_be_refined;
+ (void)initial_time;
+ (void)old_level;
- tbox::Pointer<hier::PatchHierarchy> hierarchy = patch_hierarchy;
- tbox::Pointer<geom::CartesianGridGeometry> grid_geom =
- hierarchy->getGridGeometry();
+ tbox::Pointer<hier::PatchHierarchy> hierarchy = patch_hierarchy;
+ tbox::Pointer<geom::CartesianGridGeometry> grid_geom =
+ hierarchy->getGridGeometry();
- tbox::Pointer<hier::PatchLevel> level =
- hierarchy->getPatchLevel(level_number);
+ tbox::Pointer<hier::PatchLevel> level =
+ hierarchy->getPatchLevel(level_number);
- if (allocate_data) {
- level->allocatePatchData(p_id);
- level->allocatePatchData(cell_viscosity_id);
- level->allocatePatchData(edge_viscosity_id);
- level->allocatePatchData(dp_id);
- level->allocatePatchData(p_rhs_id);
- level->allocatePatchData(p_exact_id);
- level->allocatePatchData(v_id);
- level->allocatePatchData(v_rhs_id);
- }
-
- /*
- * Initialize data in all patches in the level.
- */
- hier::PatchLevel::Iterator pi(*level);
- for (pi.initialize(*level); pi; pi++) {
-
- tbox::Pointer<hier::Patch> patch = *pi;
- if (patch.isNull()) {
- TBOX_ERROR(d_object_name
- << ": Cannot find patch. Null patch pointer.");
- }
- tbox::Pointer<pdat::CellData<double> > cell_viscosity_data =
- patch->getPatchData(cell_viscosity_id);
-
- /* At some point this needs to do the proper interpolation for
- lower levels */
- cell_viscosity_data->fill(1.0);
-
- tbox::Pointer<pdat::NodeData<double> > edge_viscosity_data =
- patch->getPatchData(edge_viscosity_id);
-
- /* At some point this needs to do the proper interpolation for
- lower levels */
- edge_viscosity_data->fill(1.0);
-
-
- tbox::Pointer<pdat::CellData<double> > dp_data =
- patch->getPatchData(dp_id);
-
- /* This is mostly so that the outer boundaries are set properly. */
- dp_data->fill(0.0);
-
- tbox::Pointer<pdat::CellData<double> > p_rhs_data =
- patch->getPatchData(p_rhs_id);
-
- p_rhs_data->fill(0.0);
-
- tbox::Pointer<pdat::SideData<double> > v_rhs_data =
- patch->getPatchData(v_rhs_id);
-
- hier::Box pbox = v_rhs_data->getBox();
-
- for(pdat::SideIterator si(pbox,0); si; si++)
- {
- pdat::SideIndex s=si();
- (*v_rhs_data)(s)=0;
- }
- for(pdat::SideIterator si(pbox,1); si; si++)
- {
- pdat::SideIndex s=si();
- (*v_rhs_data)(s)=10;
- }
- } // End patch loop.
+ if (allocate_data) {
+ level->allocatePatchData(p_id);
+ level->allocatePatchData(cell_viscosity_id);
+ level->allocatePatchData(edge_viscosity_id);
+ level->allocatePatchData(dp_id);
+ level->allocatePatchData(p_rhs_id);
+ level->allocatePatchData(p_exact_id);
+ level->allocatePatchData(v_id);
+ level->allocatePatchData(v_rhs_id);
}
+ const double inclusion_radius=0.5;
+ const double inclusion_viscosity=10;
+ const double background_viscosity=1;
+
+ /*
+ * Initialize data in all patches in the level.
+ */
+ hier::PatchLevel::Iterator pi(*level);
+ for (pi.initialize(*level); pi; pi++) {
+
+ tbox::Pointer<hier::Patch> patch = *pi;
+ if (patch.isNull()) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find patch. Null patch pointer.");
+ }
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+
+ tbox::Pointer<pdat::CellData<double> > cell_viscosity_ptr =
+ patch->getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
+
+ hier::Box visc_box = cell_viscosity.getBox();
+ for(pdat::CellIterator ci(cell_viscosity.getGhostBox()); ci; ci++)
+ {
+ pdat::CellIndex c=ci();
+ double x=geom->getXLower()[0]
+ + geom->getDx()[0]*(c[0]-visc_box.lower()[0] + 0.5);
+ double y=geom->getXLower()[1]
+ + geom->getDx()[1]*(c[1]-visc_box.lower()[1] + 0.5);
+
+ if(x*x + y*y < inclusion_radius*inclusion_radius)
+ cell_viscosity(c)=inclusion_viscosity;
+ else
+ cell_viscosity(c)=background_viscosity;
+
+ // tbox::plog << "cell "
+ // << c[0] << " "
+ // << c[1] << " "
+ // << x << " "
+ // << y << " "
+ // << geom->getXLower()[0] << " "
+ // << geom->getXLower()[1] << " "
+ // << geom->getDx()[0] << " "
+ // << geom->getDx()[1] << " "
+ // << std::boolalpha
+ // << (x*x + y*y < inclusion_radius*inclusion_radius) << " "
+ // << cell_viscosity(c) << " "
+ // << "\n";
+
+ }
+
+ tbox::Pointer<pdat::NodeData<double> > edge_viscosity_ptr =
+ patch->getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
+
+ for(pdat::NodeIterator ei(edge_viscosity.getGhostBox()); ei; ei++)
+ {
+ pdat::NodeIndex e=ei();
+ double x=geom->getXLower()[0]
+ + geom->getDx()[0]*(e[0]-visc_box.lower()[0]);
+ double y=geom->getXLower()[1]
+ + geom->getDx()[1]*(e[1]-visc_box.lower()[1]);
+ if(x*x + y*y < inclusion_radius*inclusion_radius)
+ edge_viscosity(e)=inclusion_viscosity;
+ else
+ edge_viscosity(e)=background_viscosity;
+
+ // tbox::plog << "edge "
+ // << e[0] << " "
+ // << e[1] << " "
+ // << x << " "
+ // << y << " "
+ // << geom->getXLower()[0] << " "
+ // << geom->getXLower()[1] << " "
+ // << geom->getDx()[0] << " "
+ // << geom->getDx()[1] << " "
+ // << std::boolalpha
+ // << (x*x + y*y < inclusion_radius*inclusion_radius) << " "
+ // << edge_viscosity(e) << " "
+ // << "\n";
+ }
+
+
+ // tbox::Pointer<pdat::CellData<double> > cell_viscosity_data =
+ // patch->getPatchData(cell_viscosity_id);
+
+ // /* At some point this needs to do the proper interpolation for
+ // lower levels */
+ // cell_viscosity_data->fill(1.0);
+
+ // tbox::Pointer<pdat::NodeData<double> > edge_viscosity_data =
+ // patch->getPatchData(edge_viscosity_id);
+
+ // /* At some point this needs to do the proper interpolation for
+ // lower levels */
+ // edge_viscosity_data->fill(1.0);
+
+
+ tbox::Pointer<pdat::CellData<double> > dp_data =
+ patch->getPatchData(dp_id);
+
+ /* This is mostly so that the outer boundaries are set properly. */
+ dp_data->fill(0.0);
+
+ tbox::Pointer<pdat::CellData<double> > p_rhs_data =
+ patch->getPatchData(p_rhs_id);
+
+ p_rhs_data->fill(0.0);
+
+ tbox::Pointer<pdat::SideData<double> > v_rhs_data =
+ patch->getPatchData(v_rhs_id);
+
+ hier::Box pbox = v_rhs_data->getBox();
+
+ for(pdat::SideIterator si(pbox,0); si; si++)
+ {
+ pdat::SideIndex s=si();
+ (*v_rhs_data)(s)=0;
+ }
+ for(pdat::SideIterator si(pbox,1); si; si++)
+ {
+ pdat::SideIndex s=si();
+ // (*v_rhs_data)(s)=10;
+ (*v_rhs_data)(s)=0;
+ }
+ } // End patch loop.
}
diff -r 087c7e707b4d -r c82543f37f62 StokesFACOps.h
--- a/StokesFACOps.h Sat Mar 05 11:25:10 2011 -0800
+++ b/StokesFACOps.h Sat Mar 05 11:29:15 2011 -0800
@@ -563,10 +563,9 @@ private:
pdat::NodeData<double> &edge_viscosity,
const double &theta_momentum);
- /* Compute 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). */
+ /* 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,
@@ -580,6 +579,11 @@ private:
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,
@@ -625,6 +629,43 @@ private:
center,down,right_e,center_e,dy,dx);
return result;
+ }
+
+ /* 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,
+ pdat::CellData<double> &p,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::NodeData<double> &edge_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 pdat::SideIndex &up_x,
+ const pdat::SideIndex &down_x,
+ const pdat::SideIndex ¢er_y,
+ const pdat::SideIndex &up_y,
+ const pdat::NodeIndex ¢er_e,
+ const pdat::NodeIndex &up_e,
+ const hier::Index &ip,
+ 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);
+
+ 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 087c7e707b4d -r c82543f37f62 StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Sat Mar 05 11:25:10 2011 -0800
+++ b/StokesFACOps/Update_V.C Sat Mar 05 11:29:15 2011 -0800
@@ -90,7 +90,7 @@ void SAMRAI::solv::StokesFACOps::Update_
if(!((center[axis]==pbox.lower(axis) && v(left_x)==boundary_value)
|| (center[axis]==pbox.upper(axis)+1 && v(right_x)==boundary_value)))
{
- double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
+ double d2vx_dxx, d2vx_dyy, C_vx;
/* If y==0 */
hier::Index offset(0,0);
bool set_boundary(false);
@@ -113,85 +113,57 @@ void SAMRAI::solv::StokesFACOps::Update_
dv=v(center_x+offset) - v(center_x);
}
- // const double dtau_xx_dx =
- // (v(right_x) - 2*v(center_x) + v(left_x))/(dx*dx);
-
- // const double dtau_xy_dy =
- // (v(up_x)-2*v(center_x)+v(down_x))/(dy*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);
-
- 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));
-
-
- // tbox::plog << "v "
- // << axis << " "
- // << center[0] << " "
- // << center[1] << " "
- // << dv_xx_dx << " "
- // << dtau_xx_dx << " "
- // << dv_xy_dy << " "
- // << dtau_xy_dy << " "
- // << "\n";
-
- // C_vx=-2*(1/(dx*dx) + 1/(dy*dy));
C_vx=dRm_dv(cell_viscosity,edge_viscosity,center,left,up_e,center_e,
dx,dy);
- dp_dx=(p(center)-p(left))/dx;
-
- double delta_Rx=v_rhs(center_x) - dtau_xx_dx - dtau_xy_dy + dp_dx;
+ double delta_Rx=v_rhs(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,dx,dy);
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_Rx));
+
+ {
+ 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;
+
// tbox::plog << "v " << axis << " "
- // // << v(pdat::SideIndex(center,
- // // axis,
- // // pdat::SideIndex::Lower))
- // // << " "
- // // << maxres << " "
- // // << (*v_rhs)(pdat::SideIndex(center,
- // // axis,
- // // pdat::SideIndex::Lower))
- // // << " "
- // // << &(*v_rhs)(pdat::SideIndex(center,
- // // axis,
- // // pdat::SideIndex::Lower))
- // // << " ";
- // << delta_Rx << " ";
- // // << d2vx_dxx << " "
- // // << d2vx_dyy << " "
- // // << dp_dx << " "
- // // << v(pdat::SideIndex(left,axis,
- // // pdat::SideIndex::Lower)) << " "
- // // << v(pdat::SideIndex
- // // (center,axis,
- // // pdat::SideIndex::Lower)) << " "
- // // << v(pdat::SideIndex
- // // (right,axis,
- // // pdat::SideIndex::Lower)) << " "
- // // << v(pdat::SideIndex(up,axis,
- // // pdat::SideIndex::Lower)) << " "
- // // << v(pdat::SideIndex
- // // (down,axis,
- // // pdat::SideIndex::Lower)) << " ";
-
- // // << (*p)(center) << " "
- // // << (*p)(left) << " ";
-
- // // << &(*v_rhs)(pdat::SideIndex(center,
- // // axis,
- // // pdat::SideIndex::Lower))
- // // << " "
- // // << std::boolalpha
- // // << set_boundary << " ";
+ // << center[0] << " "
+ // << center[1] << " "
+ // << v(center_x) << " "
+ // << v(left_x) << " "
+ // << v(right_x) << " "
+ // << v(up_x) << " "
+ // << v(down_x) << " "
+ // << v(up_y) << " "
+ // << v(up_y-ip) << " "
+ // << v(center_y) << " "
+ // << v(center_y-ip) << " "
+ // << v_rhs(center_x) << " "
+ // << delta_Rx << " "
+ // << p(center) << " "
+ // << p(left) << " "
+ // << edge_viscosity(up_e) << " "
+ // << edge_viscosity(center_e) << " "
+ // << cell_viscosity(center) << " "
+ // << cell_viscosity(left) << " "
+ // << dx << " "
+ // << dy << " "
+ // << dtau_xx_dx << " "
+ // << dtau_xy_dy << " "
+ // << dp_dx << " "
+ // << "\n";
+ }
v(center_x)+=delta_Rx*theta_momentum/C_vx;
diff -r 087c7e707b4d -r c82543f37f62 StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Sat Mar 05 11:25:10 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Sat Mar 05 11:29:15 2011 -0800
@@ -160,29 +160,38 @@ void SAMRAI::solv::StokesFACOps::compute
* S4. Compute residual on patches in level.
*/
+ const hier::Index ip(1,0), jp(0,1);
for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
tbox::Pointer<hier::Patch> patch = *pi;
tbox::Pointer<pdat::CellData<double> >
- p = solution.getComponentPatchData(0, *patch);
+ p_ptr = solution.getComponentPatchData(0, *patch);
+ pdat::CellData<double> &p(*p_ptr);
tbox::Pointer<pdat::SideData<double> >
- v = solution.getComponentPatchData(1, *patch);
+ v_ptr = solution.getComponentPatchData(1, *patch);
+ pdat::SideData<double> &v(*v_ptr);
tbox::Pointer<pdat::CellData<double> >
- cell_viscosity = patch->getPatchData(cell_viscosity_id);
+ cell_viscosity_ptr = patch->getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
tbox::Pointer<pdat::NodeData<double> >
- edge_viscosity = patch->getPatchData(edge_viscosity_id);
+ edge_viscosity_ptr = patch->getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
tbox::Pointer<pdat::CellData<double> >
- p_rhs = rhs.getComponentPatchData(0, *patch);
+ p_rhs_ptr = rhs.getComponentPatchData(0, *patch);
+ pdat::CellData<double> &p_rhs(*p_rhs_ptr);
tbox::Pointer<pdat::SideData<double> >
- v_rhs = rhs.getComponentPatchData(1, *patch);
+ v_rhs_ptr = rhs.getComponentPatchData(1, *patch);
+ pdat::SideData<double> &v_rhs(*v_rhs_ptr);
tbox::Pointer<pdat::CellData<double> >
- p_resid = residual.getComponentPatchData(0, *patch);
+ p_resid_ptr = residual.getComponentPatchData(0, *patch);
+ pdat::CellData<double> &p_resid(*p_resid_ptr);
tbox::Pointer<pdat::SideData<double> >
- v_resid = residual.getComponentPatchData(1, *patch);
+ v_resid_ptr = residual.getComponentPatchData(1, *patch);
+ pdat::SideData<double> &v_resid(*v_resid_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry> geom = patch->getPatchGeometry();
- double dx = *(geom->getDx());
- double dy = *(geom->getDx());
+ double dx = geom->getDx()[0];
+ double dy = geom->getDx()[1];
for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
{
@@ -200,6 +209,22 @@ void SAMRAI::solv::StokesFACOps::compute
++right[0];
--left[0];
+ const pdat::SideIndex
+ center_x(center,0,pdat::SideIndex::Lower),
+ left_x(left,0,pdat::SideIndex::Lower),
+ right_x(right,0,pdat::SideIndex::Lower),
+ down_x(down,0,pdat::SideIndex::Lower),
+ up_x(up,0,pdat::SideIndex::Lower),
+ center_y(center,1,pdat::SideIndex::Lower),
+ left_y(left,1,pdat::SideIndex::Lower),
+ right_y(right,1,pdat::SideIndex::Lower),
+ down_y(down,1,pdat::SideIndex::Lower),
+ up_y(up,1,pdat::SideIndex::Lower);
+ const pdat::NodeIndex
+ center_e(center,pdat::NodeIndex::LowerLeft),
+ up_e(up,pdat::NodeIndex::LowerLeft),
+ right_e(right,pdat::NodeIndex::LowerLeft);
+
// tbox::plog << "resid "
// << ln << " "
// << i << " "
@@ -207,155 +232,57 @@ void SAMRAI::solv::StokesFACOps::compute
/* p */
if(i!=pbox.upper(0)+1 && j!=pbox.upper(1)+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;
- (*p_resid)(center)=(*p_rhs)(center) - dvx_dx - dvy_dy;
-
+ double dvx_dx=(v(right_x) - v(center_x))/dx;
+ double dvy_dy=(v(up_y) - v(center_y))/dy;
+ p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy;
// tbox::plog << "p "
- // // << (*p)(center) << " ";
- // << (*p_resid)(center) << " ";
- // // << (*p_rhs)(center) << " "
- // // << dvx_dx << " "
- // // << dvy_dy << " "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // // pdat::SideIndex::Upper)) << " "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // // pdat::SideIndex::Lower)) << " "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- // // pdat::SideIndex::Upper)) << " "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- // // pdat::SideIndex::Lower)) << " ";
+ // << p_resid(center) << " ";
}
/* vx */
if(j!=pbox.upper(1)+1)
{
/* If x==0 */
- if((center[0]==pbox.lower(0)
- && (*v)(pdat::SideIndex(left,
- pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- ==boundary_value)
+ if((center[0]==pbox.lower(0) && v(left_x)==boundary_value)
|| (center[0]==pbox.upper(0)+1
- && (*v)(pdat::SideIndex(right,
- pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- ==boundary_value))
+ && v(right_x)==boundary_value))
+
{
- (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=0;
+ v_resid(center_x)=0;
}
else
{
- double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
- d2vx_dyy=
- ((*v)(pdat::SideIndex(up,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(down,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dy*dy);
-
- C_vx=-2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy));
-
- d2vx_dxx=((*v)(pdat::SideIndex(left,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(right,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dx*dx);
-
- dp_dx=((*p)(center)-(*p)(left))/dx;
-
- (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - (*cell_viscosity)(center)*(d2vx_dxx + d2vx_dyy) + dp_dx;
+ v_resid(center_x)=v_rhs(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,dx,dy);
+ // tbox::plog << "vx "
+ // << v_resid(center_x) << " "
+ // << v(center_x) << " "
+ // << v(right_x) << " "
+ // << v(left_x) << " ";
}
-
-
-
- // tbox::plog << "vx "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // // pdat::SideIndex::Lower))
- // // << " ";
- // << (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
- // pdat::SideIndex::Lower))
- // << " ";
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // // pdat::SideIndex::Upper))
- // // << " "
- // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // // pdat::SideIndex::Lower))
- // // << " "
- // // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // // pdat::SideIndex::Upper)))
- // // << " "
}
/* vy */
if(i!=pbox.upper(0)+1)
{
/* If y==0 */
- if((center[1]==pbox.lower(1)
- && (*v)(pdat::SideIndex(down,
- pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- ==boundary_value)
- || (center[1]==pbox.upper(1)+1
- && (*v)(pdat::SideIndex(up,
- pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- ==boundary_value))
+ if((center[1]==pbox.lower(1) && v(down_y)==boundary_value)
+ || (center[1]==pbox.upper(1)+1 && v(up_y)==boundary_value))
{
- (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=0;
+ v_resid(center_y)=0;
}
else
{
- double dp_dy, d2vy_dxx, d2vy_dyy, C_vy;
- d2vy_dxx=
- ((*v)(pdat::SideIndex(left,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(right,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dx*dx);
-
- C_vy=-2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy));
- d2vy_dyy=((*v)(pdat::SideIndex(up,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(down,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dy*dy);
-
- dp_dy=((*p)(center)-(*p)(down))/dy;
-
- (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - (*cell_viscosity)(center)*(d2vy_dxx + d2vy_dyy) + dp_dy;
+ v_resid(center_y)=v_rhs(center_y)
+ - 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,dy,dx);
}
// tbox::plog << "vy "
- // << (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::Y,
- // pdat::SideIndex::Lower))
- // << " ";
+ // << v_resid(center_y) << " ";
}
// tbox::plog << "\n";
}
More information about the CIG-COMMITS
mailing list