[cig-commits] commit: Add in node viscosity. Still assumes an isoviscous solve
Mercurial
hg at geodynamics.org
Fri Feb 25 14:17:26 PST 2011
changeset: 96:a11e70ab6421
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 23 19:29:28 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:
Add in node viscosity. Still assumes an isoviscous solve
diff -r 872d7c48020c -r a11e70ab6421 FACStokes.h
--- a/FACStokes.h Wed Feb 23 16:20:27 2011 -0800
+++ b/FACStokes.h Wed Feb 23 19:29:28 2011 -0800
@@ -168,7 +168,8 @@ namespace SAMRAI {
*
* These are initialized in the constructor and never change.
*/
- int p_id, viscosity_id, dp_id, p_exact_id, p_rhs_id, v_id, v_rhs_id;
+ int p_id, cell_viscosity_id, node_viscosity_id, dp_id, p_exact_id,
+ p_rhs_id, v_id, v_rhs_id;
//@}
diff -r 872d7c48020c -r a11e70ab6421 FACStokes/FACStokes.C
--- a/FACStokes/FACStokes.C Wed Feb 23 16:20:27 2011 -0800
+++ b/FACStokes/FACStokes.C Wed Feb 23 19:29:28 2011 -0800
@@ -70,12 +70,25 @@ namespace SAMRAI {
stencil widths */);
tbox::Pointer<pdat::CellVariable<double> >
- viscosity(new pdat::CellVariable<double>(dim,
- object_name + ":viscosity"));
- viscosity_id = vdb->registerVariableAndContext(viscosity,d_context,
- hier::IntVector(dim, 1)
- /* ghost cell width is
- 1 in case needed */);
+ cell_viscosity(new pdat::CellVariable<double>(dim,
+ object_name
+ + ":cell_viscosity"));
+ cell_viscosity_id = vdb->registerVariableAndContext(cell_viscosity,
+ d_context,
+ hier::IntVector(dim, 1)
+ /* ghost cell width is
+ 1 in case needed */);
+
+ tbox::Pointer<pdat::NodeVariable<double> >
+ node_viscosity(new pdat::NodeVariable<double>(dim,
+ object_name
+ + ":node_viscosity"));
+ node_viscosity_id = vdb->registerVariableAndContext(node_viscosity,
+ d_context,
+ hier::IntVector(dim, 1)
+ /* ghost cell width is
+ 1 in case needed */);
+
tbox::Pointer<pdat::CellVariable<double> >
dp(new pdat::CellVariable<double>(dim, object_name + ":dp"));
diff -r 872d7c48020c -r a11e70ab6421 FACStokes/initializeLevelData.C
--- a/FACStokes/initializeLevelData.C Wed Feb 23 16:20:27 2011 -0800
+++ b/FACStokes/initializeLevelData.C Wed Feb 23 19:29:28 2011 -0800
@@ -66,7 +66,8 @@ namespace SAMRAI {
if (allocate_data) {
level->allocatePatchData(p_id);
- level->allocatePatchData(viscosity_id);
+ level->allocatePatchData(cell_viscosity_id);
+ level->allocatePatchData(node_viscosity_id);
level->allocatePatchData(dp_id);
level->allocatePatchData(p_rhs_id);
level->allocatePatchData(p_exact_id);
@@ -85,12 +86,20 @@ namespace SAMRAI {
TBOX_ERROR(d_object_name
<< ": Cannot find patch. Null patch pointer.");
}
- tbox::Pointer<pdat::CellData<double> > viscosity_data =
- patch->getPatchData(viscosity_id);
+ 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 */
- viscosity_data->fill(1.0);
+ cell_viscosity_data->fill(1.0);
+
+ tbox::Pointer<pdat::NodeData<double> > node_viscosity_data =
+ patch->getPatchData(node_viscosity_id);
+
+ /* At some point this needs to do the proper interpolation for
+ lower levels */
+ node_viscosity_data->fill(1.0);
+
tbox::Pointer<pdat::CellData<double> > dp_data =
patch->getPatchData(dp_id);
diff -r 872d7c48020c -r a11e70ab6421 FACStokes/setupPlotter.C
--- a/FACStokes/setupPlotter.C Wed Feb 23 16:20:27 2011 -0800
+++ b/FACStokes/setupPlotter.C Wed Feb 23 19:29:28 2011 -0800
@@ -46,9 +46,12 @@ namespace SAMRAI {
plotter.registerPlotQuantity("Exact solution",
"SCALAR",
p_exact_id);
- plotter.registerPlotQuantity("Viscosity",
+ plotter.registerPlotQuantity("Cell Viscosity",
"SCALAR",
- viscosity_id);
+ cell_viscosity_id);
+ plotter.registerPlotQuantity("Node Viscosity",
+ "SCALAR",
+ node_viscosity_id);
plotter.registerPlotQuantity("Stokes source",
"SCALAR",
p_rhs_id);
diff -r 872d7c48020c -r a11e70ab6421 FACStokes/solveStokes.C
--- a/FACStokes/solveStokes.C Wed Feb 23 16:20:27 2011 -0800
+++ b/FACStokes/solveStokes.C Wed Feb 23 19:29:28 2011 -0800
@@ -56,12 +56,13 @@ int SAMRAI::FACStokes::solveStokes()
}
d_stokes_fac_solver.initializeSolverState
- (p_id,viscosity_id,dp_id,p_rhs_id,v_id,v_rhs_id,d_hierarchy,0,
- d_hierarchy->getFinestLevelNumber());
+ (p_id,cell_viscosity_id,node_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,viscosity_id,dp_id,
+ solver_ret = d_stokes_fac_solver.solveSystem(p_id,cell_viscosity_id,
+ node_viscosity_id,dp_id,
p_rhs_id,v_id,v_rhs_id);
/*
* Present data on the solve.
diff -r 872d7c48020c -r a11e70ab6421 StokesFACOps.h
--- a/StokesFACOps.h Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACOps.h Wed Feb 23 19:29:28 2011 -0800
@@ -296,9 +296,11 @@ public:
* flux and you would like that to be used, set flux id to the
* patch data index of that space.
*/
- void set_viscosity_dp_id(const int &viscosity, const int &dp)
+ void set_viscosity_dp_id(const int &cell_viscosity,
+ const int &node_viscosity, const int &dp)
{
- viscosity_id=viscosity;
+ cell_viscosity_id=cell_viscosity;
+ node_viscosity_id=node_viscosity;
dp_id=dp;
}
//@}
@@ -545,7 +547,8 @@ private:
double &maxres,
const double &dx,
const double &dy,
- const double &viscosity,
+ tbox::Pointer<pdat::CellData<double> > &cell_viscosity,
+ tbox::Pointer<pdat::NodeData<double> > &node_viscosity,
const double &theta_momentum);
/*!
@@ -874,7 +877,7 @@ private:
*
* @see set_viscosity_dp_id.
*/
- int viscosity_id, dp_id;
+ int cell_viscosity_id, node_viscosity_id, dp_id;
#ifdef HAVE_HYPRE
/*!
diff -r 872d7c48020c -r a11e70ab6421 StokesFACOps/StokesFACOps.C
--- a/StokesFACOps/StokesFACOps.C Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACOps/StokesFACOps.C Wed Feb 23 19:29:28 2011 -0800
@@ -97,7 +97,8 @@ namespace SAMRAI {
d_coarse_solver_tolerance(1.e-8),
d_coarse_solver_max_iterations(10),
d_residual_tolerance_during_smoothing(-1.0),
- viscosity_id(-1),
+ cell_viscosity_id(-1),
+ node_viscosity_id(-1),
dp_id(-1),
#ifdef HAVE_HYPRE
d_hypre_solver(dim,
diff -r 872d7c48020c -r a11e70ab6421 StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACOps/Update_V.C Wed Feb 23 19:29:28 2011 -0800
@@ -63,7 +63,8 @@ void SAMRAI::solv::StokesFACOps::Update_
double &maxres,
const double &dx,
const double &dy,
- const double &viscosity,
+ tbox::Pointer<pdat::CellData<double> > &cell_viscosity,
+ tbox::Pointer<pdat::NodeData<double> > &node_viscosity,
const double &theta_momentum)
{
const int off_axis=(axis==0) ? 1 : 0;
@@ -123,7 +124,9 @@ void SAMRAI::solv::StokesFACOps::Update_
pdat::SideIndex::Lower)))
/(dy*dy);
- C_vx=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
+ 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));
d2vx_dxx=((*v)(pdat::SideIndex
(left,axis,
@@ -142,7 +145,7 @@ void SAMRAI::solv::StokesFACOps::Update_
(*v_rhs)(pdat::SideIndex(center,
axis,
pdat::SideIndex::Lower))
- - viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
+ - (*cell_viscosity)(center)*(d2vx_dxx + d2vx_dyy) + dp_dx;
/* No scaling here, though there should be. */
diff -r 872d7c48020c -r a11e70ab6421 StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 23 19:29:28 2011 -0800
@@ -159,7 +159,6 @@ void SAMRAI::solv::StokesFACOps::compute
/*
* S4. Compute residual on patches in level.
*/
- double viscosity=1;
for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
tbox::Pointer<hier::Patch> patch = *pi;
@@ -167,6 +166,10 @@ void SAMRAI::solv::StokesFACOps::compute
p = solution.getComponentPatchData(0, *patch);
tbox::Pointer<pdat::SideData<double> >
v = solution.getComponentPatchData(1, *patch);
+ 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_rhs = rhs.getComponentPatchData(0, *patch);
tbox::Pointer<pdat::SideData<double> >
@@ -263,7 +266,7 @@ void SAMRAI::solv::StokesFACOps::compute
pdat::SideIndex::Lower)))
/(dy*dy);
- C_vx=-2*viscosity*(1/(dx*dx) + 1/(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))
@@ -279,7 +282,7 @@ void SAMRAI::solv::StokesFACOps::compute
pdat::SideIndex::Lower))=
(*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))
- - viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
+ - (*cell_viscosity)(center)*(d2vx_dxx + d2vx_dyy) + dp_dx;
}
@@ -332,7 +335,7 @@ void SAMRAI::solv::StokesFACOps::compute
pdat::SideIndex::Lower)))
/(dx*dx);
- C_vy=-2*viscosity*(1/(dx*dx) + 1/(dy*dy));
+ 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,
@@ -347,7 +350,7 @@ void SAMRAI::solv::StokesFACOps::compute
pdat::SideIndex::Lower))=
(*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Lower))
- - viscosity*(d2vy_dxx + d2vy_dyy) + dp_dy;
+ - (*cell_viscosity)(center)*(d2vy_dxx + d2vy_dyy) + dp_dy;
}
// tbox::plog << "vy "
// << (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::Y,
diff -r 872d7c48020c -r a11e70ab6421 StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 23 19:29:28 2011 -0800
@@ -90,7 +90,6 @@ void SAMRAI::solv::StokesFACOps::smoothE
xeqScheduleGhostFill(p_id, v_id, ln);
}
- double viscosity=1;
double theta_momentum=0.7;
double theta_continuity=1.0;
@@ -128,6 +127,10 @@ void SAMRAI::solv::StokesFACOps::smoothE
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);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -205,7 +208,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
|| (i==pbox.upper(0)+1 && j<pbox.upper(1)+1))
{
Update_V(0,j,pbox,geom,center,left,right,down,up,p,
- v,v_rhs,maxres,dx,dy,viscosity,theta_momentum);
+ v,v_rhs,maxres,dx,dy,cell_viscosity,
+ node_viscosity,theta_momentum);
}
}
}
@@ -231,6 +235,10 @@ void SAMRAI::solv::StokesFACOps::smoothE
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);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -308,7 +316,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
|| (i<pbox.upper(0)+1 && j==pbox.upper(1)+1))
{
Update_V(1,i,pbox,geom,center,down,up,left,right,p,
- v,v_rhs,maxres,dy,dx,viscosity,theta_momentum);
+ v,v_rhs,maxres,dy,dx,cell_viscosity,
+ node_viscosity,theta_momentum);
}
}
}
@@ -337,6 +346,10 @@ void SAMRAI::solv::StokesFACOps::smoothE
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);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -430,7 +443,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
maxres=std::max(maxres,std::fabs(delta_R_continuity));
(*dp)(center)=
- viscosity*delta_R_continuity*theta_continuity;
+ (*cell_viscosity)(center)*delta_R_continuity*theta_continuity;
// if(ln==2 && i==15)
@@ -480,6 +493,10 @@ void SAMRAI::solv::StokesFACOps::smoothE
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);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -569,7 +586,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
==boundary_value)))
(*v)(pdat::SideIndex(center,0,pdat::SideIndex::Lower))
-=((*dp)(center) - (*dp)(left))
- /(dx*2*viscosity*(1/(dx*dx) + 1/(dy*dy)));
+ /(dx*2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy)));
}
if(set_p[(i-gbox.lower(0))
+ (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))]
@@ -588,7 +605,7 @@ void SAMRAI::solv::StokesFACOps::smoothE
==boundary_value)))
(*v)(pdat::SideIndex(center,1,pdat::SideIndex::Lower))
-=((*dp)(center) - (*dp)(down))
- /(dy*2*viscosity*(1/(dx*dx) + 1/(dy*dy)));
+ /(dy*2*(*cell_viscosity)(center)*(1/(dx*dx) + 1/(dy*dy)));
}
}
}
diff -r 872d7c48020c -r a11e70ab6421 StokesFACSolver.h
--- a/StokesFACSolver.h Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACSolver.h Wed Feb 23 19:29:28 2011 -0800
@@ -180,7 +180,8 @@ public:
bool
solveSystem(
const int p,
- const int viscosity,
+ const int cell_viscosity,
+ const int node_viscosity,
const int dp,
const int p_rhs,
const int v,
@@ -210,7 +211,8 @@ public:
* @see solveSystem( const int, const int, tbox::Pointer< hier::PatchHierarchy >, int, int);
*/
bool
- solveSystem(const int p, const int viscosity, const int dp, const int p_rhs,
+ solveSystem(const int p, const int cell_viscosity, const int node_viscosity,
+ const int dp, const int p_rhs,
const int v, const int v_rhs);
/*!
@@ -508,7 +510,8 @@ public:
*/
void
initializeSolverState(const int p,
- const int viscosity,
+ const int cell_viscosity,
+ const int node_viscosity,
const int dp,
const int p_rhs,
const int v,
diff -r 872d7c48020c -r a11e70ab6421 StokesFACSolver/initializeSolverState.C
--- a/StokesFACSolver/initializeSolverState.C Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACSolver/initializeSolverState.C Wed Feb 23 19:29:28 2011 -0800
@@ -30,7 +30,8 @@
void SAMRAI::solv::StokesFACSolver::initializeSolverState
(const int p,
- const int viscosity,
+ const int cell_viscosity,
+ const int node_viscosity,
const int dp,
const int p_rhs,
const int v,
@@ -103,7 +104,7 @@ void SAMRAI::solv::StokesFACSolver::init
d_fac_ops.setStokesSpecifications(d_stokes_spec);
- d_fac_ops.set_viscosity_dp_id(viscosity,dp);
+ d_fac_ops.set_viscosity_dp_id(cell_viscosity,node_viscosity,dp);
createVectorWrappers(p, p_rhs, v, v_rhs);
diff -r 872d7c48020c -r a11e70ab6421 StokesFACSolver/solveSystem.C
--- a/StokesFACSolver/solveSystem.C Wed Feb 23 16:20:27 2011 -0800
+++ b/StokesFACSolver/solveSystem.C Wed Feb 23 19:29:28 2011 -0800
@@ -28,7 +28,9 @@
*************************************************************************
*/
-bool SAMRAI::solv::StokesFACSolver::solveSystem(const int p, const int viscosity,
+bool SAMRAI::solv::StokesFACSolver::solveSystem(const int p,
+ const int cell_viscosity,
+ const int node_viscosity,
const int dp, const int p_rhs,
const int v, const int v_rhs)
{
@@ -82,7 +84,8 @@ bool SAMRAI::solv::StokesFACSolver::solv
bool SAMRAI::solv::StokesFACSolver::solveSystem
(const int p,
- const int viscosity,
+ const int cell_viscosity,
+ const int node_viscosity,
const int dp,
const int p_rhs,
const int v,
@@ -116,11 +119,12 @@ bool SAMRAI::solv::StokesFACSolver::solv
<< "specified.\n");
}
#endif
- initializeSolverState(p, viscosity, dp, p_rhs, v, v_rhs,
+ initializeSolverState(p, cell_viscosity, node_viscosity, dp, p_rhs, v, v_rhs,
hierarchy, coarse_ln, fine_ln);
bool solver_rval;
- solver_rval = solveSystem(p, viscosity, dp, p_rhs, v, v_rhs);
+ solver_rval = solveSystem(p, cell_viscosity, node_viscosity,
+ dp, p_rhs, v, v_rhs);
deallocateSolverState();
More information about the CIG-COMMITS
mailing list