[cig-commits] commit: Shear corner input works with hard coded velocity boundaries. Required a lot of modifications to support non-homogeneous BCs as well as a fixes to V_Boundary_Refine. Also made max residual work even if the residual is negative.
Mercurial
hg at geodynamics.org
Fri Feb 25 14:17:13 PST 2011
changeset: 93:12d6bfe07e49
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 16 14:11:56 2011 -0800
files: FACStokes/solveStokes.C StokesFACOps.h StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/prolongErrorAndCorrect.C StokesFACOps/restrictSolution.C StokesFACOps/set_boundaries.C StokesFACOps/smoothErrorByRedBlack.C StokesFACSolver.h V_Boundary_Refine.h V_Boundary_Refine/Update_V.C V_Boundary_Refine/refine.C V_Coarsen_Patch_Strategy.h V_Refine_Patch_Strategy.C input/shear_corner.input set_V_boundary.C set_V_boundary.h
description:
Shear corner input works with hard coded velocity boundaries. Required a lot of modifications to support non-homogeneous BCs as well as a fixes to V_Boundary_Refine. Also made max residual work even if the residual is negative.
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 FACStokes/solveStokes.C
--- a/FACStokes/solveStokes.C Wed Feb 16 14:04:05 2011 -0800
+++ b/FACStokes/solveStokes.C Wed Feb 16 14:11:56 2011 -0800
@@ -52,7 +52,7 @@ int SAMRAI::FACStokes::solveStokes()
v = patch->getPatchData(v_id);
v->fill(0.0);
}
- d_stokes_fac_solver.set_boundaries(v_id,level);
+ d_stokes_fac_solver.set_boundaries(v_id,level,false);
}
d_stokes_fac_solver.initializeSolverState
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps.h
--- a/StokesFACOps.h Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps.h Wed Feb 16 14:11:56 2011 -0800
@@ -486,10 +486,19 @@ public:
void set_boundaries(const int &v_id, const int &l)
{
+ set_boundaries(v_id,l,true);
+ }
+ void set_boundaries(const int &v_id, const int &l, const bool &rhs)
+ {
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(l);
- set_boundaries(v_id,level);
+ set_boundaries(v_id,level,rhs);
}
- void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level);
+ void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level)
+ {
+ set_boundaries(v_id,level,true);
+ }
+ void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level,
+ const bool &rhs);
//@}
@@ -523,6 +532,7 @@ private:
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,
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps/Update_V.C Wed Feb 16 14:11:56 2011 -0800
@@ -51,6 +51,7 @@ void SAMRAI::solv::StokesFACOps::Update_
(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,
@@ -86,18 +87,20 @@ void SAMRAI::solv::StokesFACOps::Update_
/* If y==0 */
hier::Index offset(0,0);
bool set_boundary(false);
- if(center[axis]==pbox.lower(axis)+1)
+ if(center[axis]==pbox.lower(axis)+1
+ && !geom->getTouchesRegularBoundary(axis,0))
{
offset[axis]=-2;
set_boundary=true;
}
- else if(center[axis]==pbox.upper(axis))
+ else if(center[axis]==pbox.upper(axis)
+ && !geom->getTouchesRegularBoundary(axis,1))
{
offset[axis]=2;
set_boundary=true;
}
- double dv;
+ double dv(0);
if(set_boundary)
{
dv=(*v)(pdat::SideIndex
@@ -108,7 +111,7 @@ void SAMRAI::solv::StokesFACOps::Update_
(center,axis,
pdat::SideIndex::Lower));
}
-
+
d2vx_dyy=
((*v)(pdat::SideIndex(up,axis,
pdat::SideIndex::Lower))
@@ -142,20 +145,50 @@ void SAMRAI::solv::StokesFACOps::Update_
- viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
- tbox::plog << "v " << axis << " "
- << (*v)(pdat::SideIndex(center,
- axis,
- pdat::SideIndex::Lower))
- << " "
- << (*v_rhs)(pdat::SideIndex(center,
- axis,
- pdat::SideIndex::Lower))
- << " "
- << std::boolalpha
- << set_boundary << " ";
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,std::fabs(delta_Rx));
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,delta_Rx);
+ // 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 << " ";
(*v)(pdat::SideIndex(center,axis,
pdat::SideIndex::Lower))+=
@@ -168,9 +201,21 @@ void SAMRAI::solv::StokesFACOps::Update_
(*v)(pdat::SideIndex(center+offset,axis,
pdat::SideIndex::Lower))=
(*v)(pdat::SideIndex(center,axis,pdat::SideIndex::Lower)) + dv;
- tbox::plog << offset(0) << " "
- << offset(1) << " "
- << dv << " ";
+ // tbox::plog << "setbc "
+ // << (center+offset)(0) << " "
+ // << (center+offset)(1) << " "
+ // // << center(0) << " "
+ // // << center(1) << " "
+ // // << offset(0) << " "
+ // // << offset(1) << " "
+ // // << pbox.lower(0) << " "
+ // // << pbox.upper(0) << " "
+ // // << pbox.lower(1) << " "
+ // // << pbox.upper(1) << " "
+ // << (*v)(pdat::SideIndex(center+offset,axis,
+ // pdat::SideIndex::Lower)) << " "
+ // << (*v)(pdat::SideIndex(center,axis,pdat::SideIndex::Lower)) << " "
+ // << dv << " ";
}
}
}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 16 14:11:56 2011 -0800
@@ -117,7 +117,7 @@ void SAMRAI::solv::StokesFACOps::compute
/* S1. Fill solution ghost data. */
- set_boundaries(v_id,ln);
+ set_boundaries(v_id,ln,error_equation_indicator);
if (ln > d_ln_min) {
/* Fill from current, next coarser level and physical boundary */
xeqScheduleGhostFill(p_id, v_id, ln);
@@ -197,6 +197,10 @@ void SAMRAI::solv::StokesFACOps::compute
++right[0];
--left[0];
+ tbox::plog << "resid "
+ << ln << " "
+ << i << " "
+ << j << " ";
/* p */
if(i!=pbox.upper(0)+1 && j!=pbox.upper(1)+1)
{
@@ -211,6 +215,14 @@ void SAMRAI::solv::StokesFACOps::compute
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
pdat::SideIndex::Lower)))/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 << " ";
}
/* vx */
@@ -264,28 +276,22 @@ void SAMRAI::solv::StokesFACOps::compute
- tbox::plog << "resid "
- << ln << " "
- << i << " "
- << j << " "
- // << (*p_resid)(center) << " "
- // << (*p_rhs)(center) << " "
- // << dvx_dx << " "
- // << dvy_dy << " "
- << "vx "
- << (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
+ tbox::plog << "vx "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))
- << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Upper))
- << " "
+ << " ";
+ // << (*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)))
- << " "
- << "\n";
+ // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Upper)))
+ // << " "
}
/* vy */
@@ -335,7 +341,12 @@ void SAMRAI::solv::StokesFACOps::compute
pdat::SideIndex::Lower))
- viscosity*(d2vy_dxx + d2vy_dyy) + dp_dy;
}
+ tbox::plog << "vy "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ << " ";
}
+ tbox::plog << "\n";
}
}
@@ -359,7 +370,7 @@ void SAMRAI::solv::StokesFACOps::compute
/* We also need to set the boundaries of the rhs so that coarsening
works correctly. */
const int v_rhs_id = rhs.getComponentDescriptorIndex(1);
- set_boundaries(v_rhs_id,ln);
+ set_boundaries(v_rhs_id,ln,true);
xeqScheduleGhostFillNoCoarse(invalid_id, v_rhs_id, ln);
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps/prolongErrorAndCorrect.C
--- a/StokesFACOps/prolongErrorAndCorrect.C Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps/prolongErrorAndCorrect.C Wed Feb 16 14:11:56 2011 -0800
@@ -98,7 +98,7 @@ void SAMRAI::solv::StokesFACOps::prolong
d_side_scratch_id,
dest_ln);
- set_boundaries(s.getComponentDescriptorIndex(1),fine_level);
+ set_boundaries(s.getComponentDescriptorIndex(1),fine_level,true);
/*
* Add the refined error in the scratch space to the error currently
* residing in the destination level.
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps/restrictSolution.C
--- a/StokesFACOps/restrictSolution.C Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps/restrictSolution.C Wed Feb 16 14:11:56 2011 -0800
@@ -67,7 +67,7 @@ void SAMRAI::solv::StokesFACOps::restric
xeqScheduleURestriction(p_dst,p_src,v_dst,v_src,dest_ln);
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(dest_ln);
- set_boundaries(v_dst,level);
+ set_boundaries(v_dst,level,false);
// v_refine_patch_strategy.setHomogeneousBc(false);
p_refine_patch_strategy.setTargetDataId(d.getComponentDescriptorIndex(0));
v_refine_patch_strategy.setTargetDataId(d.getComponentDescriptorIndex(1));
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps/set_boundaries.C
--- a/StokesFACOps/set_boundaries.C Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps/set_boundaries.C Wed Feb 16 14:11:56 2011 -0800
@@ -46,11 +46,11 @@
/* Set the physical boundaries for the velocity. */
void SAMRAI::solv::StokesFACOps::set_boundaries
-(const int &v_id, tbox::Pointer<hier::PatchLevel> &level)
+(const int &v_id, tbox::Pointer<hier::PatchLevel> &level, const bool &rhs)
{
for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
{
tbox::Pointer<hier::Patch> patch = *pi;
- set_V_boundary(*patch,v_id);
+ set_V_boundary(*patch,v_id,rhs);
}
}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 16 14:11:56 2011 -0800
@@ -78,8 +78,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
correct. */
p_refine_patch_strategy.setTargetDataId(p_id);
v_refine_patch_strategy.setTargetDataId(v_id);
+ set_boundaries(v_id,level,true);
xeqScheduleGhostFillNoCoarse(p_rhs_id,v_rhs_id,ln);
- set_boundaries(v_id,level);
if (ln > d_ln_min) {
/*
@@ -237,21 +237,16 @@ void SAMRAI::solv::StokesFACOps::smoothE
++right[0];
--left[0];
- tbox::plog << "smooth "
- << ln << " "
- << sweep << " "
- << rb << " "
- << i << " "
- << j << " "
- << std::boolalpha
- << set_p[(i-gbox.lower(0))
- + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1))] << " "
- << (i-gbox.lower(0))
- + (gbox.upper(0)-gbox.lower(0)+1)*(j-gbox.lower(1)) << " ";
- // << pbox.lower(0) << " "
- // << pbox.upper(0) << " "
- // << pbox.lower(1) << " "
- // << pbox.upper(1) << " ";
+ // tbox::plog << "smooth "
+ // << ln << " "
+ // << sweep << " "
+ // << rb << " "
+ // << i << " "
+ // << j << " ";
+ // // << pbox.lower(0) << " "
+ // // << pbox.upper(0) << " "
+ // // << pbox.lower(1) << " "
+ // // << pbox.upper(1) << " ";
/* Update p */
if(set_p[(i-gbox.lower(0))
@@ -272,43 +267,49 @@ void SAMRAI::solv::StokesFACOps::smoothE
(*p_rhs)(center) - dvx_dx - dvy_dy;
/* No scaling here, though there should be. */
- maxres=std::max(maxres,delta_R_continuity);
+ maxres=std::max(maxres,std::fabs(delta_R_continuity));
(*p)(center)+=
viscosity*delta_R_continuity*theta_continuity;
- tbox::plog << "p "
- << (*p)(center) << " "
- << maxres << " ";
- // << (*p_rhs)(center) << " "
- // << dvx_dx << " "
- // << dvy_dy << " "
- // << delta_R_continuity << " "
- // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // pdat::SideIndex::Upper)) << " "
- // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // pdat::SideIndex::Lower)) << " "
- // << dx << " ";
+ // tbox::plog << "p "
+ // // << (*p)(center) << " "
+ // // << maxres << " "
+ // // << (*p_rhs)(center) << " "
+ // // << dvx_dx << " "
+ // // << dvy_dy << " "
+ // << delta_R_continuity << " ";
+ // // << (*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)) << " ";
+ // // << dx << " "
+ // // << dy << " ";
}
/* 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))
{
- Update_V(0,j,pbox,center,left,right,down,up,p,v,v_rhs,
- maxres,dx,dy,viscosity,theta_momentum);
+ Update_V(0,j,pbox,geom,center,left,right,down,up,p,
+ v,v_rhs,maxres,dx,dy,viscosity,theta_momentum);
}
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))
{
- Update_V(1,i,pbox,center,down,up,left,right,p,v,v_rhs,
- maxres,dy,dx,viscosity,theta_momentum);
+ Update_V(1,i,pbox,geom,center,down,up,left,right,p,
+ v,v_rhs,maxres,dy,dx,viscosity,theta_momentum);
}
- tbox::plog << "\n";
+ // tbox::plog << "\n";
}
}
}
+ set_boundaries(v_id,level,true);
}
// if (residual_tolerance >= 0.0) {
/*
@@ -327,8 +328,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
converged=(tmp==1);
if (d_enable_logging)
tbox::plog
- << d_object_name << "\n"
- << " RBGS sweep #" << sweep << " : " << maxres << "\n";
+ // << d_object_name << "\n"
+ << "Smooth " << ln << " " << sweep << " : " << maxres << "\n";
// }
}
}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 StokesFACSolver.h
--- a/StokesFACSolver.h Wed Feb 16 14:04:05 2011 -0800
+++ b/StokesFACSolver.h Wed Feb 16 14:11:56 2011 -0800
@@ -558,9 +558,11 @@ public:
double
getResidualNorm() const;
- void set_boundaries(const int &v_id, tbox::Pointer<hier::PatchLevel> &level)
+ void set_boundaries(const int &v_id,
+ tbox::Pointer<hier::PatchLevel> &level,
+ const bool &homogeneous)
{
- d_fac_ops.set_boundaries(v_id,level);
+ d_fac_ops.set_boundaries(v_id,level,homogeneous);
}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 V_Boundary_Refine.h
--- a/V_Boundary_Refine.h Wed Feb 16 14:04:05 2011 -0800
+++ b/V_Boundary_Refine.h Wed Feb 16 14:11:56 2011 -0800
@@ -141,6 +141,7 @@ private:
const hier::Index &ip, const hier::Index &jp,
int &i,
int &j,
+ const int &i_max,
const int &j_max,
tbox::Pointer<pdat::SideData<double> > &v,
tbox::Pointer<pdat::SideData<double> > &v_fine) const;
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 V_Boundary_Refine/Update_V.C
--- a/V_Boundary_Refine/Update_V.C Wed Feb 16 14:04:05 2011 -0800
+++ b/V_Boundary_Refine/Update_V.C Wed Feb 16 14:11:56 2011 -0800
@@ -33,6 +33,7 @@ void SAMRAI::geom::V_Boundary_Refine::Up
const hier::Index &ip, const hier::Index &jp,
int &i,
int &j,
+ const int &i_max,
const int &j_max,
tbox::Pointer<SAMRAI::pdat::SideData<double> > &v,
tbox::Pointer<SAMRAI::pdat::SideData<double> > &v_fine) const
@@ -41,6 +42,12 @@ void SAMRAI::geom::V_Boundary_Refine::Up
center.coarsen(hier::Index(2,2));
const int off_axis= (axis==0) ? 1 : 0;
+ // tbox::plog << "Boundary Update V "
+ // << axis << " "
+ // << fine[0] << " "
+ // << fine[1] << " "
+ // << center[0] << " "
+ // << center[1] << " ";
/* Set the derivative for the normal direction */
if(boundary_direction==axis)
{
@@ -61,36 +68,47 @@ void SAMRAI::geom::V_Boundary_Refine::Up
dv_fine_plus,dv_fine_minus);
hier::Index offset(ip);
- if(!boundary_positive)
+
+ if(boundary_positive)
+ {
+ offset[axis]=2;
+ }
+ else
{
offset[axis]=-2;
+ dv_fine_minus=-dv_fine_minus;
+ dv_fine_plus=-dv_fine_plus;
}
-
- (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_minus;
- (*v_fine)(fine+jp)=(*v_fine)(fine-offset+jp) + dv_fine_plus;
-
- // tbox::plog << "Update V "
- // << axis << " "
- // << fine[0] << " "
- // << fine[1] << " "
- // << center[0] << " "
- // << center[1] << " "
+
+ if(j%2==0)
+ {
+ (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_minus/2;
+ (*v_fine)(fine+jp)=(*v_fine)(fine-offset+jp) + dv_fine_plus/2;
+ /* Since we update two points on j at once, we increment j
+ again. This is ok, since the box in the 'i' direction is
+ defined to be only one cell wide */
+ ++j;
+ }
+ else
+ {
+ (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_plus/2;
+ }
+ // tbox::plog << offset[0] << " "
+ // << offset[1] << " "
// << (*v_fine)(fine) << " "
// << (*v_fine)(fine+jp) << " "
+ // << (*v_fine)(fine-offset) << " "
+ // << (*v_fine)(fine-offset+jp) << " "
// << (*v)(center) << " "
// << (*v)(center+jp) << " "
// << (*v)(center-jp) << " "
- // << (&(*v)(center-jp)) << " "
- // << v_plus << " "
- // << v_minus << " "
- // << dv_plus << " "
- // << dv_minus << " "
+ // // << (&(*v)(center-jp)) << " "
+ // // << v_plus << " "
+ // // << v_minus << " "
+ // << dv_fine_plus << " "
+ // << dv_fine_minus << " "
// << "\n";
- /* Since we update two points on j at once, we increment j again.
- This is ok, since the box in the 'i' direction is defined to be
- only one cell wide */
- ++j;
}
/* Set the value for the tangential direction */
else
@@ -99,17 +117,42 @@ void SAMRAI::geom::V_Boundary_Refine::Up
v_center=
quad_offset_interpolate((*v)(center-jp),(*v)(center),(*v)(center+jp));
- (*v_fine)(fine)=v_center;
- if(j<j_max)
+ if(i%2==0)
+ {
+ (*v_fine)(fine)=v_center;
+
+ // tbox::plog << (*v_fine)(fine) << " ";
+ // // << (*v_fine)(fine+jp) << " "
+ // // << (*v_fine)(fine-offset) << " "
+ // // << (*v_fine)(fine-offset+jp) << " "
+ // // << (*v)(center) << " "
+ // // << (*v)(center+jp) << " "
+ // // << (*v)(center-jp) << " "
+ // // << dv_fine_plus << " "
+ // // << dv_fine_minus << " "
+
+ if(i<i_max)
+ {
+ v_plus=quad_offset_interpolate((*v)(center+ip-jp),(*v)(center+ip),
+ (*v)(center+ip+jp));
+ (*v_fine)(fine+ip)=(v_center+v_plus)/2;
+
+ // tbox::plog << (*v_fine)(fine+ip) << " ";
+
+ /* Since we update two points on 'i' at once, we increment 'i' again.
+ This is ok, since the box in the 'j' direction is defined to be
+ only one cell wide */
+ ++i;
+ }
+ }
+ else
{
v_plus=quad_offset_interpolate((*v)(center+ip-jp),(*v)(center+ip),
(*v)(center+ip+jp));
- (*v_fine)(fine+ip)=(v_center+v_plus)/2;
+ (*v_fine)(fine)=(v_center+v_plus)/2;
+ // tbox::plog << (*v_fine)(fine) << " ";
}
- /* Since we update two points on 'i' at once, we increment 'i' again.
- This is ok, since the box in the 'j' direction is defined to be
- only one cell wide */
- ++i;
+ // tbox::plog << "\n";
}
}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 V_Boundary_Refine/refine.C
--- a/V_Boundary_Refine/refine.C Wed Feb 16 14:04:05 2011 -0800
+++ b/V_Boundary_Refine/refine.C Wed Feb 16 14:11:56 2011 -0800
@@ -127,30 +127,30 @@ void SAMRAI::geom::V_Boundary_Refine::re
}
}
- tbox::plog << "VBR "
- << fine.getPatchLevelNumber() << " "
- << axis << " "
- << boundary_direction << " "
- << std::boolalpha
- << boundary_positive << " "
- << coarse_box.lower(0) << " "
- << coarse_box.upper(0) << " "
- << coarse_box.lower(1) << " "
- << coarse_box.upper(1) << " "
- << fine_box.lower(0) << " "
- << fine_box.upper(0) << " "
- << fine_box.lower(1) << " "
- << fine_box.upper(1) << " "
+ // tbox::plog << "VBR "
+ // << fine.getPatchLevelNumber() << " "
+ // << axis << " "
+ // << boundary_direction << " "
+ // << std::boolalpha
+ // << boundary_positive << " "
+ // << coarse_box.lower(0) << " "
+ // << coarse_box.upper(0) << " "
+ // << coarse_box.lower(1) << " "
+ // << coarse_box.upper(1) << " "
+ // << fine_box.lower(0) << " "
+ // << fine_box.upper(0) << " "
+ // << fine_box.lower(1) << " "
+ // << fine_box.upper(1) << " "
- << overlap_box.lower(0) << " "
- << overlap_box.upper(0) << " "
- << overlap_box.lower(1) << " "
- << overlap_box.upper(1) << " "
- << i_min << " "
- << i_max << " "
- << j_min << " "
- << j_max << " "
- << "\n";
+ // << overlap_box.lower(0) << " "
+ // << overlap_box.upper(0) << " "
+ // << overlap_box.lower(1) << " "
+ // << overlap_box.upper(1) << " "
+ // << i_min << " "
+ // << i_max << " "
+ // << j_min << " "
+ // << j_max << " "
+ // << "\n";
for(int j=j_min; j<=j_max; ++j)
for(int i=i_min; i<=i_max; ++i)
@@ -161,12 +161,12 @@ void SAMRAI::geom::V_Boundary_Refine::re
if(axis==0)
{
Update_V(axis,boundary_direction,boundary_positive,fine,
- ip,jp,i,j,j_max,v,v_fine);
+ ip,jp,i,j,i_max,j_max,v,v_fine);
}
else if(axis==1)
{
Update_V(axis,boundary_direction,boundary_positive,fine,
- jp,ip,j,i,i_max,v,v_fine);
+ jp,ip,j,i,j_max,i_max,v,v_fine);
}
}
}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 V_Coarsen_Patch_Strategy.h
--- a/V_Coarsen_Patch_Strategy.h Wed Feb 16 14:04:05 2011 -0800
+++ b/V_Coarsen_Patch_Strategy.h Wed Feb 16 14:11:56 2011 -0800
@@ -69,7 +69,7 @@ public:
const hier::Box& coarse_box,
const hier::IntVector& ratio)
{
- set_V_boundary(fine,v_id);
+ set_V_boundary(fine,v_id,true);
}
virtual void
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 V_Refine_Patch_Strategy.C
--- a/V_Refine_Patch_Strategy.C Wed Feb 16 14:04:05 2011 -0800
+++ b/V_Refine_Patch_Strategy.C Wed Feb 16 14:11:56 2011 -0800
@@ -27,7 +27,7 @@ SAMRAI::solv::V_Refine_Patch_Strategy::p
// << pbox.upper(1) << " "
// << "\n";
- set_V_boundary(coarse,v_id);
+ set_V_boundary(coarse,v_id,true);
// for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 input/shear_corner.input
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/input/shear_corner.input Wed Feb 16 14:11:56 2011 -0800
@@ -0,0 +1,148 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Input file for example FAC Stokes solver
+ *
+ ************************************************************************/
+// This is the input file for the 2D FAC example
+// demonstrating changes in boundary conditions.
+//
+// Note that using constant refine prolongation
+// lead to slower convergence.
+
+Main {
+ dim = 2
+ // Base name for output files.
+ base_name = "shear_corner"
+ // Whether to log all nodes in a parallel run.
+ log_all_nodes = TRUE
+ // Visualization writers to write files for.
+ vis_writer = "Vizamrai", "VisIt"
+}
+
+FACStokes {
+ // The FACStokes class is the "user class" in this example.
+ // It owns the solver and contains the code to set up the solver.
+ // The inputs for FACStokes is simply the inputs for the individual
+ // parts owned by the FACStokes class.
+ fac_solver {
+ // This is the input for the cell-centered Stokes FAC solver
+ // class in the SAMRAI library.
+ enable_logging = TRUE // Bool flag to switch logging on/off
+ max_cycles = 100 // Max number of FAC cycles to use
+ residual_tol = 1e-8 // Residual tolerance to solve for
+ num_pre_sweeps = 5 // Number of presmoothing sweeps to use
+ num_post_sweeps = 5 // Number of postsmoothing sweeps to use
+ coarse_solver_max_iterations = 25
+ coarse_solver_tolerance = 1e-10
+ p_prolongation_method = "P_REFINE" // Type of refinement
+ // used in prolongation.
+ // Suggested values are
+ // "LINEAR_REFINE"
+ // "CONSTANT_REFINE"
+ v_prolongation_method = "V_REFINE" // Type of refinement
+ // used in prolongation.
+ // Suggested values are
+ // "LINEAR_REFINE"
+ // "CONSTANT_REFINE"
+ use_smg = TRUE // Whether to use HYPRE's SMG instead of PFMG.
+ }
+ bc_coefs {
+ // These are the boundary condition specifications. The number
+ // after "boundary_" is the location index of the boundary.
+ // The inputs are arrays of strings where the first string
+ // indicates the type of values you want to set. "slope" means
+ // boundary slope, "value" means boundary value, and "coefficients"
+ // mean the raw Robin boundary condition coefficients.
+ // The remaining strings are converted into numbers as
+ // appropriate for what boundary condition you specified with
+ // the first string. Other boundary conditions are possible.
+ // see the solv_RobinBcCoefStrategy class.
+ // Examples:
+ // boundary_0 = "slope", "0"
+ // boundary_1 = "coefficients", "0", "0"
+ // boundary_2 = "value", "0"
+ // boundary_3 = "value", "0"
+ boundary_0 = "value", "0"
+ boundary_1 = "value", "0"
+ boundary_2 = "value", "0"
+ boundary_3 = "value", "0"
+ }
+}
+
+CartesianGridGeometry {
+ // Specify lower/upper corners of the computational domain and a
+ // set of non-overlapping boxes defining domain interior. If union
+ // of boxes is not a parallelpiped, lower/upper corner data corresponds
+ // to min/max corner indices over all boxes given.
+ // x_lo -- (double array) lower corner of computational domain [REQD]
+ // x_up -- (double array) upper corner of computational domain [REQD]
+ // domain_boxes -- (box array) set of boxes that define interior of
+ // physical domain. [REQD]
+ // periodic_dimension -- (int array) coordinate directions in which
+ // domain is periodic. Zero indicates not
+ // periodic, non-zero value indicates periodicity.
+ // [0]
+ domain_boxes = [(0,0), (3,3)]
+ x_lo = 0, 0
+ x_up = 1, 1
+}
+
+StandardTagAndInitialize {
+ tagging_method = "REFINE_BOXES"
+ RefineBoxes {
+ level_0 = [(0,0),(3,3)]
+ // level_1 = [(0,0),(7,7)]
+ // level_2 = [(0,0),(15,15)]
+ level_1 = [(2,2),(3,3)]
+ level_2 = [(4,4),(5,5)]
+ // level_3 = [(8,8),(20,20)]
+ //etc.
+ }
+}
+
+PatchHierarchy {
+ // Information used to create patches in AMR hierarchy.
+ // max_levels -- (int) max number of mesh levels in hierarchy [REQD]
+ //
+ // For most of the following parameters, the number of precribed data
+ // values need not match the number of levels in the hierarchy
+ // (determined by max_levels). If more values are given than number
+ // of levels, extraneous values will be ignored. If less are give, then
+ // values that correspond to individual levels will apply to those
+ // levels. Missing values will be taken from those for the finest
+ // level specified.
+ //
+ // ratio_to_coarser {
+ // level_1 -- (int array) ratio between index spaces on
+ // level 1 to level 0 [REQD]
+ // level_2 -- (int array) ratio between index spaces on
+ // level 2 to level 1 [REQD]
+ // etc....
+ // }
+ // largest_patch_size {
+ // level_0 -- (int array) largest patch allowed on level 0.
+ // [REQD]
+ // level_1 -- (int array) " " " " level 1
+ // [level 0 entry]
+ // etc....
+ // }
+ max_levels = 3
+ ratio_to_coarser {
+ level_1 = 2, 2
+ level_2 = 2, 2
+ level_3 = 2, 2
+ level_4 = 2, 2
+ }
+ largest_patch_size {
+ level_0 = 32, 32
+ // level_0 = 4, 4
+ // all finer levels will use same values as level_0...
+ }
+}
+
+GriddingAlgorithm {
+}
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 set_V_boundary.C
--- a/set_V_boundary.C Wed Feb 16 14:04:05 2011 -0800
+++ b/set_V_boundary.C Wed Feb 16 14:11:56 2011 -0800
@@ -7,7 +7,8 @@
using namespace SAMRAI;
-void set_V_boundary(const SAMRAI::hier::Patch& patch, const int &v_id)
+void set_V_boundary(const SAMRAI::hier::Patch& patch, const int &v_id,
+ const bool &rhs)
{
tbox::Pointer<pdat::SideData<double> >
v = patch.getPatchData(v_id);
@@ -29,6 +30,12 @@ void set_V_boundary(const SAMRAI::hier::
/* vx */
if(j<=gbox.upper(1))
{
+ // tbox::plog << "V boundary "
+ // << i << " "
+ // << j << " "
+ // << pbox.lower(0) << " "
+ // << pbox.upper(0) << " ";
+
/* Set a sentinel value */
if((i<pbox.lower(0) && geom->getTouchesRegularBoundary(0,0))
|| (i>pbox.upper(0)+1 && geom->getTouchesRegularBoundary(0,1)))
@@ -36,6 +43,13 @@ void set_V_boundary(const SAMRAI::hier::
(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))=
boundary_value;
+ }
+ /* Set vx=-1 on the right boundary */
+ else if(i==pbox.upper(0)+1 && geom->getTouchesRegularBoundary(0,1)
+ && !rhs)
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=-1;
}
/* Set the value so the derivative=0 */
else if(j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
@@ -52,10 +66,14 @@ void set_V_boundary(const SAMRAI::hier::
(*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
pdat::SideIndex::Lower));
}
+ // tbox::plog << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower)) << " "
+ // << "\n";
}
/* vy */
if(i<=gbox.upper(0))
{
+ /* Set a sentinel value */
if((j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
|| (j>pbox.upper(1)+1 && geom->getTouchesRegularBoundary(1,1)))
{
@@ -63,6 +81,14 @@ void set_V_boundary(const SAMRAI::hier::
pdat::SideIndex::Lower))=
boundary_value;
}
+ /* Set vy=1 on the top boundary */
+ else if(j==pbox.upper(1)+1 && geom->getTouchesRegularBoundary(1,1)
+ && !rhs)
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=1;
+ }
+ /* Set the value so the derivative=0 */
else if(i<pbox.lower(0) && geom->getTouchesRegularBoundary(0,0))
{
(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
diff -r 4d2ba42b4e79 -r 12d6bfe07e49 set_V_boundary.h
--- a/set_V_boundary.h Wed Feb 16 14:04:05 2011 -0800
+++ b/set_V_boundary.h Wed Feb 16 14:11:56 2011 -0800
@@ -3,6 +3,7 @@
#include "SAMRAI/hier/Patch.h"
-void set_V_boundary(const SAMRAI::hier::Patch& patch, const int &id);
+void set_V_boundary(const SAMRAI::hier::Patch& patch, const int &id,
+ const bool &rhs);
#endif
More information about the CIG-COMMITS
mailing list