[cig-commits] commit: Make smooth look at the residuals and stop early if needed.
Mercurial
hg at geodynamics.org
Fri Feb 25 14:15:05 PST 2011
changeset: 53:fdc7a6ec07dc
user: Walter Landry <wlandry at caltech.edu>
date: Thu Jan 13 12:54:39 2011 -0800
files: StokesFACOps/smoothErrorByRedBlack.C
description:
Make smooth look at the residuals and stop early if needed.
diff -r 277940b2cb1c -r fdc7a6ec07dc StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Thu Jan 13 12:12:27 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Thu Jan 13 12:54:39 2011 -0800
@@ -87,261 +87,272 @@ void SAMRAI::solv::StokesFACOps::smoothE
* Smooth the number of sweeps specified or until
* the convergence is satisfactory.
*/
- int isweep(0);
- double red_maxres, blk_maxres, maxres = 0;
- red_maxres = blk_maxres = residual_tolerance + 1;
+ double maxres;
/*
- * Instead of checking residual convergence globally,
- * we check the not_converged flag. This avoids possible
- * round-off errors affecting different processes differently,
- * leading to disagreement on whether to continue smoothing.
+ * Instead of checking residual convergence globally, we check the
+ * converged flag. This avoids possible round-off errors affecting
+ * different processes differently, leading to disagreement on
+ * whether to continue smoothing.
*/
bool converged = false;
- for (; isweep < num_sweeps*(1<<(d_ln_max-ln)) && !converged; ++isweep) {
- red_maxres = blk_maxres = 0;
+ for (int sweep=0; sweep < num_sweeps*(1<<(d_ln_max-ln)) && !converged; ++sweep)
+ {
+ maxres=0;
+ for(int rb=0;rb<2;++rb)
+ {
+ // Need to sync
+ xeqScheduleGhostFillNoCoarse(p_id,v_id,ln);
+ for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
+ tbox::Pointer<hier::Patch> patch = *pi;
- for(int rb=0;rb<2;++rb)
- {
- // Need to sync
- xeqScheduleGhostFillNoCoarse(p_id,v_id,ln);
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++) {
- 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> >
- 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);
+ hier::Box pbox=patch->getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch->getPatchGeometry();
+ double dx = *(geom->getDx());
+ double dy = *(geom->getDx());
- hier::Box pbox=patch->getBox();
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = patch->getPatchGeometry();
- double dx = *(geom->getDx());
- double dy = *(geom->getDx());
+ for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
+ {
+ /* Do the red-black skip */
+ int i_min=pbox.lower(0) + (abs(pbox.lower(0) + j + rb))%2;
+ for(int i=i_min; i<=pbox.upper(0)+1; i+=2)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
- for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
- {
- /* Do the red-black skip */
- int i_min=pbox.lower(0) + (abs(pbox.lower(0) + j + rb))%2;
- for(int i=i_min; i<=pbox.upper(0)+1; i+=2)
- {
- pdat::CellIndex center(tbox::Dimension(2));
- center[0]=i;
- center[1]=j;
+ pdat::CellIndex up(center), down(center), right(center),
+ left(center);
- pdat::CellIndex up(center), down(center), right(center),
- left(center);
+ ++up[1];
+ --down[1];
+ ++right[0];
+ --left[0];
- ++up[1];
- --down[1];
- ++right[0];
- --left[0];
+ /* Update 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;
- /* Update 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;
+ double delta_R_continuity=
+ (*p_rhs)(center) - dvx_dx - dvy_dy;
- double delta_R_continuity=
- (*p_rhs)(center) - dvx_dx - dvy_dy;
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,delta_R_continuity);
- (*p)(center)+=
- viscosity*delta_R_continuity*theta_continuity;
- }
+ (*p)(center)+=
+ viscosity*delta_R_continuity*theta_continuity;
+ }
- /* Update vx */
- if(j!=pbox.upper(1)+1)
- {
- /* If x==0 */
- if((center[0]==pbox.lower(0)
- && geom->getTouchesRegularBoundary(0,0))
- || (center[0]==pbox.upper(0)+1
- && geom->getTouchesRegularBoundary(0,1)))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=0;
- }
- else
- {
- double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
- /* If y==0 */
- if(center[1]==pbox.lower(1)
- && geom->getTouchesRegularBoundary(1,0))
- {
- d2vx_dyy=
- ((*v)(pdat::SideIndex(up,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex(center,
- pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dy*dy);
- C_vx=-viscosity*(2/(dx*dx) + 1/(dy*dy));
- }
- /* If y==max_y */
- else if(center[1]==pbox.upper(1)
- && geom->getTouchesRegularBoundary(1,1))
- {
- d2vx_dyy=
- (-(*v)(pdat::SideIndex(center,
+ /* Update vx */
+ if(j!=pbox.upper(1)+1)
+ {
+ /* If x==0 */
+ if((center[0]==pbox.lower(0)
+ && geom->getTouchesRegularBoundary(0,0))
+ || (center[0]==pbox.upper(0)+1
+ && geom->getTouchesRegularBoundary(0,1)))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=0;
+ }
+ else
+ {
+ double dp_dx, d2vx_dxx, d2vx_dyy, C_vx;
+ /* If y==0 */
+ if(center[1]==pbox.lower(1)
+ && geom->getTouchesRegularBoundary(1,0))
+ {
+ d2vx_dyy=
+ ((*v)(pdat::SideIndex(up,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ - (*v)(pdat::SideIndex
+ (center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))
+ /(dy*dy);
+ C_vx=-viscosity*(2/(dx*dx) + 1/(dy*dy));
+ }
+ /* If y==max_y */
+ else if(center[1]==pbox.upper(1)
+ && geom->getTouchesRegularBoundary(1,1))
+ {
+ d2vx_dyy=
+ (-(*v)(pdat::SideIndex(center,
+ pdat::SideIndex::X,
+ pdat::SideIndex::Lower))
+ + (*v)(pdat::SideIndex
+ (down,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)))
+ /(dy*dy);
+ C_vx=-viscosity*(2/(dx*dx) + 1/(dy*dy));
+ }
+ else
+ {
+ 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*viscosity*(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;
+
+ double delta_Rx=
+ (*v_rhs)(pdat::SideIndex(center,
pdat::SideIndex::X,
pdat::SideIndex::Lower))
- + (*v)(pdat::SideIndex(down,
- pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- /(dy*dy);
- C_vx=-viscosity*(2/(dx*dx) + 1/(dy*dy));
- }
- else
- {
- 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);
+ - viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
- C_vx=-2*viscosity*(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);
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,delta_Rx);
- dp_dx=((*p)(center)-(*p)(left))/dx;
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))+=
+ delta_Rx*theta_momentum/C_vx;
+ }
+ }
+
+ /* Update vy */
+ if(i!=pbox.upper(0)+1)
+ {
+ /* If y==0 */
+ if((center[1]==pbox.lower(1)
+ && geom->getTouchesRegularBoundary(1,0))
+ || (center[1]==pbox.upper(1)+1
+ && geom->getTouchesRegularBoundary(1,1)))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=0;
+ }
+ else
+ {
+ double dp_dy, d2vy_dxx, d2vy_dyy, C_vy;
+ /* If x==0 */
+ if(center[0]==pbox.lower(0)
+ && geom->getTouchesRegularBoundary(0,0))
+ {
+ d2vy_dxx=
+ ((*v)(pdat::SideIndex(right,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ - (*v)(pdat::SideIndex
+ (center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))
+ /(dx*dx);
+ C_vy=-viscosity*(1/(dx*dx) + 2/(dy*dy));
+ }
+ /* If x==max_x */
+ else if(center[0]==pbox.upper(0)
+ && geom->getTouchesRegularBoundary(0,1))
+ {
+ d2vy_dxx=
+ ((*v)(pdat::SideIndex(left,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ - (*v)(pdat::SideIndex
+ (center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower)))
+ /(dx*dx);
+ C_vy=-viscosity*(1/(dx*dx) + 2/(dy*dy));
+ }
+ else
+ {
+ 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*viscosity*(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;
- double delta_Rx=
- (*v_rhs)(pdat::SideIndex(center,
- pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- - viscosity*(d2vx_dxx + d2vx_dyy) + dp_dx;
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))+=
- delta_Rx*theta_momentum/C_vx;
- }
- }
+ double delta_Ry=
+ (*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))
+ - viscosity*(d2vy_dxx + d2vy_dyy) + dp_dy;
- /* Update vy */
- if(i!=pbox.upper(0)+1)
- {
- /* If y==0 */
- if((center[1]==pbox.lower(1)
- && geom->getTouchesRegularBoundary(1,0))
- || (center[1]==pbox.upper(1)+1
- && geom->getTouchesRegularBoundary(1,1)))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=0;
- }
- else
- {
- double dp_dy, d2vy_dxx, d2vy_dyy, C_vy;
- /* If x==0 */
- if(center[0]==pbox.lower(0)
- && geom->getTouchesRegularBoundary(0,0))
- {
- d2vy_dxx=
- ((*v)(pdat::SideIndex(right,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex(center,
- pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dx*dx);
- C_vy=-viscosity*(1/(dx*dx) + 2/(dy*dy));
- }
- /* If x==max_x */
- else if(center[0]==pbox.upper(0)
- && geom->getTouchesRegularBoundary(0,1))
- {
- d2vy_dxx=
- ((*v)(pdat::SideIndex(left,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - (*v)(pdat::SideIndex(center,
- pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- /(dx*dx);
- C_vy=-viscosity*(1/(dx*dx) + 2/(dy*dy));
- }
- else
- {
- 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);
+ /* No scaling here, though there should be. */
+ maxres=std::max(maxres,delta_Ry);
- C_vy=-2*viscosity*(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;
-
- double delta_Ry=
- (*v_rhs)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- - viscosity*(d2vy_dxx + d2vy_dyy) + dp_dy;
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))+=
- delta_Ry*theta_momentum/C_vy;
- }
- }
- }
- }
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))+=
+ delta_Ry*theta_momentum/C_vy;
+ }
+ }
+ }
+ }
+ }
}
+ if (residual_tolerance >= 0.0) {
+ /*
+ * Check for early end of sweeps due to convergence
+ * only if it is numerically possible (user gave a
+ * non negative value for residual tolerance).
+ */
+ converged = maxres < residual_tolerance;
+ const tbox::SAMRAI_MPI&
+ mpi(d_hierarchy->getDomainMappedBoxLevel().getMPI());
+ int tmp= converged ? 1 : 0;
+ if (mpi.getSize() > 1)
+ {
+ mpi.AllReduce(&tmp, 1, MPI_MIN);
+ }
+ converged=(tmp==1);
+ if (d_enable_logging)
+ tbox::plog
+ << d_object_name << "\n"
+ << " RBGS sweep #" << sweep << " : " << maxres << "\n";
}
- }
-
- // if (residual_tolerance >= 0.0) {
- // /*
- // * Check for early end of sweeps due to convergence
- // * only if it is numerically possible (user gave a
- // * non negative value for residual tolerance).
- // */
- // maxres = tbox::MathUtilities<double>::Max(red_maxres, blk_maxres);
- // not_converged = maxres > residual_tolerance;
- // const tbox::SAMRAI_MPI& mpi(d_hierarchy->getDomainMappedBoxLevel().getMPI());
- // if (mpi.getSize() > 1) {
- // mpi.AllReduce(¬_converged, 1, MPI_MAX);
- // }
- // }
-
- // if (d_enable_logging) tbox::plog
- // << d_object_name << " RBGS smoothing maxres = " << maxres << "\n"
- // << " after " << isweep << " sweeps.\n";
-
+ }
}
More information about the CIG-COMMITS
mailing list