[cig-commits] commit: Variable viscosity setup with constant viscosity seems to work, though not as well as the one set up just for constant viscosity
Mercurial
hg at geodynamics.org
Fri Feb 25 14:17:30 PST 2011
changeset: 97:b09ac3138b25
user: Walter Landry <wlandry at caltech.edu>
date: Fri Feb 25 11:47:03 2011 -0800
files: FACStokes.h FACStokes/FACStokes.C FACStokes/initializeLevelData.C FACStokes/setupPlotter.C FACStokes/solveStokes.C StokesFACOps.h StokesFACOps/StokesFACOps.C StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/smoothErrorByRedBlack.C StokesFACSolver.h StokesFACSolver/initializeSolverState.C StokesFACSolver/solveSystem.C
description:
Variable viscosity setup with constant viscosity seems to work, though not as well as the one set up just for constant viscosity
diff -r a11e70ab6421 -r b09ac3138b25 FACStokes.h
--- a/FACStokes.h Wed Feb 23 19:29:28 2011 -0800
+++ b/FACStokes.h Fri Feb 25 11:47:03 2011 -0800
@@ -168,7 +168,7 @@ namespace SAMRAI {
*
* These are initialized in the constructor and never change.
*/
- int p_id, cell_viscosity_id, node_viscosity_id, dp_id, p_exact_id,
+ int p_id, cell_viscosity_id, edge_viscosity_id, dp_id, p_exact_id,
p_rhs_id, v_id, v_rhs_id;
//@}
diff -r a11e70ab6421 -r b09ac3138b25 FACStokes/FACStokes.C
--- a/FACStokes/FACStokes.C Wed Feb 23 19:29:28 2011 -0800
+++ b/FACStokes/FACStokes.C Fri Feb 25 11:47:03 2011 -0800
@@ -80,10 +80,10 @@ namespace SAMRAI {
1 in case needed */);
tbox::Pointer<pdat::NodeVariable<double> >
- node_viscosity(new pdat::NodeVariable<double>(dim,
+ edge_viscosity(new pdat::NodeVariable<double>(dim,
object_name
- + ":node_viscosity"));
- node_viscosity_id = vdb->registerVariableAndContext(node_viscosity,
+ + ":edge_viscosity"));
+ edge_viscosity_id = vdb->registerVariableAndContext(edge_viscosity,
d_context,
hier::IntVector(dim, 1)
/* ghost cell width is
diff -r a11e70ab6421 -r b09ac3138b25 FACStokes/initializeLevelData.C
--- a/FACStokes/initializeLevelData.C Wed Feb 23 19:29:28 2011 -0800
+++ b/FACStokes/initializeLevelData.C Fri Feb 25 11:47:03 2011 -0800
@@ -67,7 +67,7 @@ namespace SAMRAI {
if (allocate_data) {
level->allocatePatchData(p_id);
level->allocatePatchData(cell_viscosity_id);
- level->allocatePatchData(node_viscosity_id);
+ level->allocatePatchData(edge_viscosity_id);
level->allocatePatchData(dp_id);
level->allocatePatchData(p_rhs_id);
level->allocatePatchData(p_exact_id);
@@ -93,12 +93,12 @@ namespace SAMRAI {
lower levels */
cell_viscosity_data->fill(1.0);
- tbox::Pointer<pdat::NodeData<double> > node_viscosity_data =
- patch->getPatchData(node_viscosity_id);
+ 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 */
- node_viscosity_data->fill(1.0);
+ edge_viscosity_data->fill(1.0);
tbox::Pointer<pdat::CellData<double> > dp_data =
diff -r a11e70ab6421 -r b09ac3138b25 FACStokes/setupPlotter.C
--- a/FACStokes/setupPlotter.C Wed Feb 23 19:29:28 2011 -0800
+++ b/FACStokes/setupPlotter.C Fri Feb 25 11:47:03 2011 -0800
@@ -49,9 +49,6 @@ namespace SAMRAI {
plotter.registerPlotQuantity("Cell Viscosity",
"SCALAR",
cell_viscosity_id);
- plotter.registerPlotQuantity("Node Viscosity",
- "SCALAR",
- node_viscosity_id);
plotter.registerPlotQuantity("Stokes source",
"SCALAR",
p_rhs_id);
diff -r a11e70ab6421 -r b09ac3138b25 FACStokes/solveStokes.C
--- a/FACStokes/solveStokes.C Wed Feb 23 19:29:28 2011 -0800
+++ b/FACStokes/solveStokes.C Fri Feb 25 11:47:03 2011 -0800
@@ -56,13 +56,13 @@ int SAMRAI::FACStokes::solveStokes()
}
d_stokes_fac_solver.initializeSolverState
- (p_id,cell_viscosity_id,node_viscosity_id,dp_id,p_rhs_id,v_id,v_rhs_id,
+ (p_id,cell_viscosity_id,edge_viscosity_id,dp_id,p_rhs_id,v_id,v_rhs_id,
d_hierarchy,0,d_hierarchy->getFinestLevelNumber());
tbox::plog << "solving..." << std::endl;
int solver_ret;
solver_ret = d_stokes_fac_solver.solveSystem(p_id,cell_viscosity_id,
- node_viscosity_id,dp_id,
+ edge_viscosity_id,dp_id,
p_rhs_id,v_id,v_rhs_id);
/*
* Present data on the solve.
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACOps.h
--- a/StokesFACOps.h Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACOps.h Fri Feb 25 11:47:03 2011 -0800
@@ -297,10 +297,10 @@ public:
* patch data index of that space.
*/
void set_viscosity_dp_id(const int &cell_viscosity,
- const int &node_viscosity, const int &dp)
+ const int &edge_viscosity, const int &dp)
{
cell_viscosity_id=cell_viscosity;
- node_viscosity_id=node_viscosity;
+ edge_viscosity_id=edge_viscosity;
dp_id=dp;
}
//@}
@@ -533,23 +533,43 @@ private:
int num_sweeps,
double residual_tolerance = -1.0);
- void Update_V(const int &axis, const int j,
- const hier::Box &pbox,
- tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ void Update_V
+ (const int &axis,
+ const int j,
+ const hier::Box &pbox,
+ tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ const pdat::CellIndex ¢er,
+ const pdat::CellIndex &left,
+ const pdat::CellIndex &right,
+ const pdat::CellIndex &down,
+ const pdat::CellIndex &up,
+ pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::SideData<double> &v_rhs,
+ double &maxres,
+ const double &dx,
+ const double &dy,
+ pdat::CellData<double> &cell_viscosity,
+ 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). */
+
+ double dRm_dv(pdat::CellData<double> &cell_viscosity,
+ pdat::NodeData<double> &edge_viscosity,
const pdat::CellIndex ¢er,
const pdat::CellIndex &left,
- const pdat::CellIndex &right,
- const pdat::CellIndex &down,
- const pdat::CellIndex &up,
- tbox::Pointer<pdat::CellData<double> > &p,
- tbox::Pointer<pdat::SideData<double> > &v,
- tbox::Pointer<pdat::SideData<double> > &v_rhs,
- double &maxres,
+ const pdat::NodeIndex &up_e,
+ const pdat::NodeIndex ¢er_e,
const double &dx,
- const double &dy,
- tbox::Pointer<pdat::CellData<double> > &cell_viscosity,
- tbox::Pointer<pdat::NodeData<double> > &node_viscosity,
- const double &theta_momentum);
+ const double &dy)
+ {
+ return -2*(cell_viscosity(center) + cell_viscosity(left))/(dx*dx)
+ - (edge_viscosity(up_e) + edge_viscosity(center_e))/(dy*dy);
+ }
/*!
* @brief Solve the coarsest level using HYPRE
@@ -877,7 +897,7 @@ private:
*
* @see set_viscosity_dp_id.
*/
- int cell_viscosity_id, node_viscosity_id, dp_id;
+ int cell_viscosity_id, edge_viscosity_id, dp_id;
#ifdef HAVE_HYPRE
/*!
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACOps/StokesFACOps.C
--- a/StokesFACOps/StokesFACOps.C Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACOps/StokesFACOps.C Fri Feb 25 11:47:03 2011 -0800
@@ -98,7 +98,7 @@ namespace SAMRAI {
d_coarse_solver_max_iterations(10),
d_residual_tolerance_during_smoothing(-1.0),
cell_viscosity_id(-1),
- node_viscosity_id(-1),
+ edge_viscosity_id(-1),
dp_id(-1),
#ifdef HAVE_HYPRE
d_hypre_solver(dim,
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACOps/Update_V.C Fri Feb 25 11:47:03 2011 -0800
@@ -57,32 +57,38 @@ void SAMRAI::solv::StokesFACOps::Update_
const pdat::CellIndex &right,
const pdat::CellIndex &down,
const pdat::CellIndex &up,
- tbox::Pointer<pdat::CellData<double> > &p,
- tbox::Pointer<pdat::SideData<double> > &v,
- tbox::Pointer<pdat::SideData<double> > &v_rhs,
+ pdat::CellData<double> &p,
+ pdat::SideData<double> &v,
+ pdat::SideData<double> &v_rhs,
double &maxres,
const double &dx,
const double &dy,
- tbox::Pointer<pdat::CellData<double> > &cell_viscosity,
- tbox::Pointer<pdat::NodeData<double> > &node_viscosity,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::NodeData<double> &edge_viscosity,
const double &theta_momentum)
{
const int off_axis=(axis==0) ? 1 : 0;
+ hier::Index ip(0,0);
+ ip[axis]=1;
+
+ const pdat::SideIndex
+ center_x(center,axis,pdat::SideIndex::Lower),
+ left_x(left,axis,pdat::SideIndex::Lower),
+ right_x(right,axis,pdat::SideIndex::Lower),
+ down_x(down,axis,pdat::SideIndex::Lower),
+ up_x(up,axis,pdat::SideIndex::Lower),
+ center_y(center,off_axis,pdat::SideIndex::Lower),
+ up_y(up,off_axis,pdat::SideIndex::Lower);
+ const pdat::NodeIndex
+ center_e(center,pdat::NodeIndex::LowerLeft),
+ up_e(up,pdat::NodeIndex::LowerLeft);
+
/* Update vx */
if(j<pbox.upper(off_axis)+1)
{
/* If at the 'x' boundaries, leave vx as is */
- if(!((center[axis]==pbox.lower(axis)
- && (*v)(pdat::SideIndex(left,
- axis,
- pdat::SideIndex::Lower))
- ==boundary_value)
- || (center[axis]==pbox.upper(axis)+1
- && (*v)(pdat::SideIndex
- (right,
- axis,
- pdat::SideIndex::Lower))
- ==boundary_value)))
+ 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;
/* If y==0 */
@@ -104,55 +110,49 @@ void SAMRAI::solv::StokesFACOps::Update_
double dv(0);
if(set_boundary)
{
- dv=(*v)(pdat::SideIndex
- (center+offset,
- axis,
- pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex
- (center,axis,
- pdat::SideIndex::Lower));
+ dv=v(center_x+offset) - v(center_x);
}
- d2vx_dyy=
- ((*v)(pdat::SideIndex(up,axis,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex
- (center,axis,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex
- (down,axis,
- pdat::SideIndex::Lower)))
- /(dy*dy);
+ // const double dv_xx_dx =
+ // (v(right_x) - 2*v(center_x) + v(left_x))/(dx*dx);
- C_vx=-2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy));
- // C_vx=-2*((*cell_viscosity)(center) + (*cell_viscosity)(left))/(dx*dx)
- // - ((*node_viscosity)(up) + (*node_viscosity)(center))/(dx*dx)/(dy*dy));
+ // const double dv_xy_dy =
+ // (v(up_x)-2*v(center_x)+v(down_x))/(dy*dy);
- d2vx_dxx=((*v)(pdat::SideIndex
- (left,axis,
- pdat::SideIndex::Lower))
- - 2*(*v)(pdat::SideIndex
- (center,axis,
- pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex
- (right,axis,
- pdat::SideIndex::Lower)))
- /(dx*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);
- dp_dx=((*p)(center)-(*p)(left))/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)(pdat::SideIndex(center,
- axis,
- pdat::SideIndex::Lower))
- - (*cell_viscosity)(center)*(d2vx_dxx + d2vx_dyy) + dp_dx;
-
+ double delta_Rx=v_rhs(center_x) - dtau_xx_dx - dtau_xy_dy + dp_dx;
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_Rx));
// tbox::plog << "v " << axis << " "
- // // << (*v)(pdat::SideIndex(center,
+ // // << v(pdat::SideIndex(center,
// // axis,
// // pdat::SideIndex::Lower))
// // << " "
@@ -169,17 +169,17 @@ void SAMRAI::solv::StokesFACOps::Update_
// // << d2vx_dxx << " "
// // << d2vx_dyy << " "
// // << dp_dx << " "
- // // << (*v)(pdat::SideIndex(left,axis,
+ // // << v(pdat::SideIndex(left,axis,
// // pdat::SideIndex::Lower)) << " "
- // // << (*v)(pdat::SideIndex
+ // // << v(pdat::SideIndex
// // (center,axis,
// // pdat::SideIndex::Lower)) << " "
- // // << (*v)(pdat::SideIndex
+ // // << v(pdat::SideIndex
// // (right,axis,
// // pdat::SideIndex::Lower)) << " "
- // // << (*v)(pdat::SideIndex(up,axis,
+ // // << v(pdat::SideIndex(up,axis,
// // pdat::SideIndex::Lower)) << " "
- // // << (*v)(pdat::SideIndex
+ // // << v(pdat::SideIndex
// // (down,axis,
// // pdat::SideIndex::Lower)) << " ";
@@ -193,17 +193,14 @@ void SAMRAI::solv::StokesFACOps::Update_
// // << std::boolalpha
// // << set_boundary << " ";
- (*v)(pdat::SideIndex(center,axis,
- pdat::SideIndex::Lower))+=
- delta_Rx*theta_momentum/C_vx;
+ v(center_x)+=delta_Rx*theta_momentum/C_vx;
/* Set the boundary elements so that the
derivative is zero. */
if(set_boundary)
{
- (*v)(pdat::SideIndex(center+offset,axis,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center,axis,pdat::SideIndex::Lower)) + dv;
+ v(center_x+offset)=v(center_x) + dv;
+
// tbox::plog << "setbc "
// << (center+offset)(0) << " "
// << (center+offset)(1) << " "
@@ -215,9 +212,9 @@ void SAMRAI::solv::StokesFACOps::Update_
// // << pbox.upper(0) << " "
// // << pbox.lower(1) << " "
// // << pbox.upper(1) << " "
- // << (*v)(pdat::SideIndex(center+offset,axis,
+ // << v(pdat::SideIndex(center+offset,axis,
// pdat::SideIndex::Lower)) << " "
- // << (*v)(pdat::SideIndex(center,axis,pdat::SideIndex::Lower)) << " "
+ // << v(pdat::SideIndex(center,axis,pdat::SideIndex::Lower)) << " "
// << dv << " ";
}
}
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Fri Feb 25 11:47:03 2011 -0800
@@ -169,7 +169,7 @@ void SAMRAI::solv::StokesFACOps::compute
tbox::Pointer<pdat::CellData<double> >
cell_viscosity = patch->getPatchData(cell_viscosity_id);
tbox::Pointer<pdat::NodeData<double> >
- node_viscosity = patch->getPatchData(node_viscosity_id);
+ edge_viscosity = patch->getPatchData(edge_viscosity_id);
tbox::Pointer<pdat::CellData<double> >
p_rhs = rhs.getComponentPatchData(0, *patch);
tbox::Pointer<pdat::SideData<double> >
@@ -222,18 +222,18 @@ void SAMRAI::solv::StokesFACOps::compute
// 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) << " ";
+ // // << (*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)) << " ";
}
/* vx */
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Fri Feb 25 11:47:03 2011 -0800
@@ -93,6 +93,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
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
* the convergence is satisfactory.
@@ -119,18 +121,29 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
tbox::Pointer<hier::Patch> patch = *pi;
- tbox::Pointer<pdat::CellData<double> >
- p = patch->getPatchData(p_id);
- tbox::Pointer<pdat::CellData<double> >
- p_rhs = patch->getPatchData(p_rhs_id);
- tbox::Pointer<pdat::SideData<double> >
- v = patch->getPatchData(v_id);
- tbox::Pointer<pdat::SideData<double> >
- v_rhs = patch->getPatchData(v_rhs_id);
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity = patch->getPatchData(cell_viscosity_id);
- tbox::Pointer<pdat::NodeData<double> >
- node_viscosity = patch->getPatchData(node_viscosity_id);
+ tbox::Pointer<pdat::CellData<double> > p_ptr =
+ patch->getPatchData(p_id);
+ pdat::CellData<double> &p(*p_ptr);
+ tbox::Pointer<pdat::CellData<double> > dp_ptr =
+ patch->getPatchData(dp_id);
+ pdat::CellData<double> &dp(*dp_ptr);
+ tbox::Pointer<pdat::CellData<double> > p_rhs_ptr =
+ patch->getPatchData(p_rhs_id);
+ pdat::CellData<double> &p_rhs(*p_rhs_ptr);
+
+ tbox::Pointer<pdat::SideData<double> > v_ptr =
+ patch->getPatchData(v_id);
+ pdat::SideData<double> &v(*v_ptr);
+ tbox::Pointer<pdat::SideData<double> > v_rhs_ptr =
+ patch->getPatchData(v_rhs_id);
+ pdat::SideData<double> &v_rhs(*v_rhs_ptr);
+
+ tbox::Pointer<pdat::CellData<double> > cell_visc_ptr
+ = patch->getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_visc_ptr);
+ tbox::Pointer<pdat::NodeData<double> > edge_visc_ptr
+ = patch->getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -142,7 +155,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
should set the pressure or just let it be. This is
needed at coarse/fine boundaries where the pressure
is fixed. */
- hier::Box gbox=p->getGhostBox();
+ hier::Box gbox=p.getGhostBox();
std::vector<bool> set_p(gbox.size(),true);
const tbox::Array<hier::BoundaryBox >&edges
@@ -209,7 +222,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
Update_V(0,j,pbox,geom,center,left,right,down,up,p,
v,v_rhs,maxres,dx,dy,cell_viscosity,
- node_viscosity,theta_momentum);
+ edge_viscosity,theta_momentum);
}
}
}
@@ -227,18 +240,29 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
tbox::Pointer<hier::Patch> patch = *pi;
- tbox::Pointer<pdat::CellData<double> >
- p = patch->getPatchData(p_id);
- tbox::Pointer<pdat::CellData<double> >
- p_rhs = patch->getPatchData(p_rhs_id);
- tbox::Pointer<pdat::SideData<double> >
- v = patch->getPatchData(v_id);
- tbox::Pointer<pdat::SideData<double> >
- v_rhs = patch->getPatchData(v_rhs_id);
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity = patch->getPatchData(cell_viscosity_id);
- tbox::Pointer<pdat::NodeData<double> >
- node_viscosity = patch->getPatchData(node_viscosity_id);
+ tbox::Pointer<pdat::CellData<double> > p_ptr =
+ patch->getPatchData(p_id);
+ pdat::CellData<double> &p(*p_ptr);
+ tbox::Pointer<pdat::CellData<double> > dp_ptr =
+ patch->getPatchData(dp_id);
+ pdat::CellData<double> &dp(*dp_ptr);
+ tbox::Pointer<pdat::CellData<double> > p_rhs_ptr =
+ patch->getPatchData(p_rhs_id);
+ pdat::CellData<double> &p_rhs(*p_rhs_ptr);
+
+ tbox::Pointer<pdat::SideData<double> > v_ptr =
+ patch->getPatchData(v_id);
+ pdat::SideData<double> &v(*v_ptr);
+ tbox::Pointer<pdat::SideData<double> > v_rhs_ptr =
+ patch->getPatchData(v_rhs_id);
+ pdat::SideData<double> &v_rhs(*v_rhs_ptr);
+
+ tbox::Pointer<pdat::CellData<double> > cell_visc_ptr
+ = patch->getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_visc_ptr);
+ tbox::Pointer<pdat::NodeData<double> > edge_visc_ptr
+ = patch->getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -250,7 +274,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
should set the pressure or just let it be. This is
needed at coarse/fine boundaries where the pressure
is fixed. */
- hier::Box gbox=p->getGhostBox();
+ hier::Box gbox=p.getGhostBox();
std::vector<bool> set_p(gbox.size(),true);
const tbox::Array<hier::BoundaryBox >&edges
@@ -317,7 +341,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
Update_V(1,i,pbox,geom,center,down,up,left,right,p,
v,v_rhs,maxres,dy,dx,cell_viscosity,
- node_viscosity,theta_momentum);
+ edge_viscosity,theta_momentum);
}
}
}
@@ -336,20 +360,29 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
tbox::Pointer<hier::Patch> patch = *pi;
- tbox::Pointer<pdat::CellData<double> >
- p = patch->getPatchData(p_id);
- tbox::Pointer<pdat::CellData<double> >
- dp = patch->getPatchData(dp_id);
- tbox::Pointer<pdat::CellData<double> >
- p_rhs = patch->getPatchData(p_rhs_id);
- tbox::Pointer<pdat::SideData<double> >
- v = patch->getPatchData(v_id);
- tbox::Pointer<pdat::SideData<double> >
- v_rhs = patch->getPatchData(v_rhs_id);
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity = patch->getPatchData(cell_viscosity_id);
- tbox::Pointer<pdat::NodeData<double> >
- node_viscosity = patch->getPatchData(node_viscosity_id);
+ tbox::Pointer<pdat::CellData<double> > p_ptr =
+ patch->getPatchData(p_id);
+ pdat::CellData<double> &p(*p_ptr);
+ tbox::Pointer<pdat::CellData<double> > dp_ptr =
+ patch->getPatchData(dp_id);
+ pdat::CellData<double> &dp(*dp_ptr);
+ tbox::Pointer<pdat::CellData<double> > p_rhs_ptr =
+ patch->getPatchData(p_rhs_id);
+ pdat::CellData<double> &p_rhs(*p_rhs_ptr);
+
+ tbox::Pointer<pdat::SideData<double> > v_ptr =
+ patch->getPatchData(v_id);
+ pdat::SideData<double> &v(*v_ptr);
+ tbox::Pointer<pdat::SideData<double> > v_rhs_ptr =
+ patch->getPatchData(v_rhs_id);
+ pdat::SideData<double> &v_rhs(*v_rhs_ptr);
+
+ tbox::Pointer<pdat::CellData<double> > cell_visc_ptr
+ = patch->getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_visc_ptr);
+ tbox::Pointer<pdat::NodeData<double> > edge_visc_ptr
+ = patch->getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -361,7 +394,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
should set the pressure or just let it be. This is
needed at coarse/fine boundaries where the pressure
is fixed. */
- hier::Box gbox=p->getGhostBox();
+ hier::Box gbox=p.getGhostBox();
std::vector<bool> set_p(gbox.size(),true);
const tbox::Array<hier::BoundaryBox >&edges
@@ -421,52 +454,88 @@ void SAMRAI::solv::StokesFACOps::smoothE
++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);
+
/* 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,
+ (v(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Upper))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ - v(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower)))/dx;
double dvy_dy=
- ((*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ (v(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Upper))
- - (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ - v(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Lower)))/dy;
double delta_R_continuity=
- (*p_rhs)(center) - dvx_dx - dvy_dy;
+ p_rhs(center) - dvx_dx - dvy_dy;
/* No scaling here, though there should be. */
maxres=std::max(maxres,std::fabs(delta_R_continuity));
- (*dp)(center)=
- (*cell_viscosity)(center)*delta_R_continuity*theta_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);
- // if(ln==2 && i==15)
+ 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;
+
+
+
+ // if(ln==2)
// tbox::plog << "smooth p "
// << i << " "
// << j << " "
- // << (*dp)(center) << " "
- // << delta_R_continuity << " "
- // << dvx_dx << " "
- // << dvy_dy << " "
- // << (*p_rhs)(center) << " "
- // << (*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)) << " "
+ // << dRc_dp << " "
+ // // << dp(center) << " "
+ // // << delta_R_continuity << " "
+ // // << dvx_dx << " "
+ // // << dvy_dy << " "
+ // // << p_rhs(center) << " "
+ // // << 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)) << " "
// << "\n";
- (*p)(center)+=(*dp)(center);
+ p(center)+=dp(center);
}
}
}
@@ -483,20 +552,29 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
tbox::Pointer<hier::Patch> patch = *pi;
- tbox::Pointer<pdat::CellData<double> >
- p = patch->getPatchData(p_id);
- tbox::Pointer<pdat::CellData<double> >
- dp = patch->getPatchData(dp_id);
- tbox::Pointer<pdat::CellData<double> >
- p_rhs = patch->getPatchData(p_rhs_id);
- tbox::Pointer<pdat::SideData<double> >
- v = patch->getPatchData(v_id);
- tbox::Pointer<pdat::SideData<double> >
- v_rhs = patch->getPatchData(v_rhs_id);
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity = patch->getPatchData(cell_viscosity_id);
- tbox::Pointer<pdat::NodeData<double> >
- node_viscosity = patch->getPatchData(node_viscosity_id);
+ tbox::Pointer<pdat::CellData<double> > p_ptr =
+ patch->getPatchData(p_id);
+ pdat::CellData<double> &p(*p_ptr);
+ tbox::Pointer<pdat::CellData<double> > dp_ptr =
+ patch->getPatchData(dp_id);
+ pdat::CellData<double> &dp(*dp_ptr);
+ tbox::Pointer<pdat::CellData<double> > p_rhs_ptr =
+ patch->getPatchData(p_rhs_id);
+ pdat::CellData<double> &p_rhs(*p_rhs_ptr);
+
+ tbox::Pointer<pdat::SideData<double> > v_ptr =
+ patch->getPatchData(v_id);
+ pdat::SideData<double> &v(*v_ptr);
+ tbox::Pointer<pdat::SideData<double> > v_rhs_ptr =
+ patch->getPatchData(v_rhs_id);
+ pdat::SideData<double> &v_rhs(*v_rhs_ptr);
+
+ tbox::Pointer<pdat::CellData<double> > cell_visc_ptr
+ = patch->getPatchData(cell_viscosity_id);
+ pdat::CellData<double> &cell_viscosity(*cell_visc_ptr);
+ tbox::Pointer<pdat::NodeData<double> > edge_visc_ptr
+ = patch->getPatchData(edge_viscosity_id);
+ pdat::NodeData<double> &edge_viscosity(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -508,7 +586,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
should set the pressure or just let it be. This is
needed at coarse/fine boundaries where the pressure
is fixed. */
- hier::Box gbox=p->getGhostBox();
+ hier::Box gbox=p.getGhostBox();
std::vector<bool> set_p(gbox.size(),true);
const tbox::Array<hier::BoundaryBox >&edges
@@ -568,44 +646,64 @@ void SAMRAI::solv::StokesFACOps::smoothE
++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),
+ center_y(center,1,pdat::SideIndex::Lower),
+ up_y(up,1,pdat::SideIndex::Lower),
+ down_y(down,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);
+
/* Update v */
if(set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
|| (i==pbox.upper(0)+1 && j<pbox.upper(1)+1))
{
if(!((center[0]==pbox.lower(0)
- && (*v)(pdat::SideIndex(left,
- 0,
- pdat::SideIndex::Lower))
- ==boundary_value)
+ && v(left_x)==boundary_value)
|| (center[0]==pbox.upper(0)+1
- && (*v)(pdat::SideIndex
- (right,
- 0,
- pdat::SideIndex::Lower))
- ==boundary_value)))
- (*v)(pdat::SideIndex(center,0,pdat::SideIndex::Lower))
- -=((*dp)(center) - (*dp)(left))
- /(dx*2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy)));
+ && v(right_x)==boundary_value)))
+
+ // v(pdat::SideIndex(center,0,pdat::SideIndex::Lower))
+ // -=(dp(center) - dp(left))
+ // /(dx*3*(1/(dx*dx) + 1/(dy*dy)));
+
+ // tbox::plog << "dRm_dv "
+ // << i << " "
+ // << j << " "
+ // << -3*(1/(dx*dx) + 1/(dy*dy)) << " "
+ // << dRm_dv(cell_viscosity,edge_viscosity,center,
+ // left,up_e,center_e,dx,dy) << " "
+ // << (dp(center) - dp(left))
+ // /(dx*2*(1/(dx*dx) + 1/(dy*dy))) << " "
+ // << (dp(center) - dp(left))
+ // /(dx*dRm_dv(cell_viscosity,edge_viscosity,center,
+ // left,up_e,center_e,dx,dy)) << " "
+ // << "\n";
+
+ v(center_x)+=(dp(center) - dp(left))
+ /(dx*dRm_dv(cell_viscosity,edge_viscosity,center,
+ left,up_e,center_e,dx,dy));
}
if(set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
|| (i<pbox.upper(0)+1 && j==pbox.upper(1)+1))
{
if(!((center[1]==pbox.lower(1)
- && (*v)(pdat::SideIndex(down,
- 1,
- pdat::SideIndex::Lower))
- ==boundary_value)
+ && v(down_y)==boundary_value)
|| (center[1]==pbox.upper(1)+1
- && (*v)(pdat::SideIndex
- (up,
- 1,
- pdat::SideIndex::Lower))
- ==boundary_value)))
- (*v)(pdat::SideIndex(center,1,pdat::SideIndex::Lower))
- -=((*dp)(center) - (*dp)(down))
- /(dy*2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy)));
+ && v(up_y)==boundary_value)))
+
+ // v(pdat::SideIndex(center,1,pdat::SideIndex::Lower))
+ // -=(dp(center) - dp(down))
+ // /(dy*3*(1/(dx*dx) + 1/(dy*dy)));
+ v(center_y)+=(dp(center) - dp(down))
+ /(dy*dRm_dv(cell_viscosity,edge_viscosity,center,
+ down,right_e,center_e,dy,dx));
}
}
}
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACSolver.h
--- a/StokesFACSolver.h Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACSolver.h Fri Feb 25 11:47:03 2011 -0800
@@ -181,7 +181,7 @@ public:
solveSystem(
const int p,
const int cell_viscosity,
- const int node_viscosity,
+ const int edge_viscosity,
const int dp,
const int p_rhs,
const int v,
@@ -211,7 +211,7 @@ public:
* @see solveSystem( const int, const int, tbox::Pointer< hier::PatchHierarchy >, int, int);
*/
bool
- solveSystem(const int p, const int cell_viscosity, const int node_viscosity,
+ solveSystem(const int p, const int cell_viscosity, const int edge_viscosity,
const int dp, const int p_rhs,
const int v, const int v_rhs);
@@ -511,7 +511,7 @@ public:
void
initializeSolverState(const int p,
const int cell_viscosity,
- const int node_viscosity,
+ const int edge_viscosity,
const int dp,
const int p_rhs,
const int v,
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACSolver/initializeSolverState.C
--- a/StokesFACSolver/initializeSolverState.C Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACSolver/initializeSolverState.C Fri Feb 25 11:47:03 2011 -0800
@@ -31,7 +31,7 @@ void SAMRAI::solv::StokesFACSolver::init
void SAMRAI::solv::StokesFACSolver::initializeSolverState
(const int p,
const int cell_viscosity,
- const int node_viscosity,
+ const int edge_viscosity,
const int dp,
const int p_rhs,
const int v,
@@ -104,7 +104,7 @@ void SAMRAI::solv::StokesFACSolver::init
d_fac_ops.setStokesSpecifications(d_stokes_spec);
- d_fac_ops.set_viscosity_dp_id(cell_viscosity,node_viscosity,dp);
+ d_fac_ops.set_viscosity_dp_id(cell_viscosity,edge_viscosity,dp);
createVectorWrappers(p, p_rhs, v, v_rhs);
diff -r a11e70ab6421 -r b09ac3138b25 StokesFACSolver/solveSystem.C
--- a/StokesFACSolver/solveSystem.C Wed Feb 23 19:29:28 2011 -0800
+++ b/StokesFACSolver/solveSystem.C Fri Feb 25 11:47:03 2011 -0800
@@ -30,7 +30,7 @@
bool SAMRAI::solv::StokesFACSolver::solveSystem(const int p,
const int cell_viscosity,
- const int node_viscosity,
+ const int edge_viscosity,
const int dp, const int p_rhs,
const int v, const int v_rhs)
{
@@ -85,7 +85,7 @@ bool SAMRAI::solv::StokesFACSolver::solv
bool SAMRAI::solv::StokesFACSolver::solveSystem
(const int p,
const int cell_viscosity,
- const int node_viscosity,
+ const int edge_viscosity,
const int dp,
const int p_rhs,
const int v,
@@ -119,11 +119,11 @@ bool SAMRAI::solv::StokesFACSolver::solv
<< "specified.\n");
}
#endif
- initializeSolverState(p, cell_viscosity, node_viscosity, dp, p_rhs, v, v_rhs,
+ initializeSolverState(p, cell_viscosity, edge_viscosity, dp, p_rhs, v, v_rhs,
hierarchy, coarse_ln, fine_ln);
bool solver_rval;
- solver_rval = solveSystem(p, cell_viscosity, node_viscosity,
+ solver_rval = solveSystem(p, cell_viscosity, edge_viscosity,
dp, p_rhs, v, v_rhs);
deallocateSolverState();
More information about the CIG-COMMITS
mailing list