[cig-commits] commit: Make it work in 3D
Mercurial
hg at geodynamics.org
Thu Jun 7 13:35:49 PDT 2012
changeset: 265:8e656c4c1517
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Tue May 01 15:15:15 2012 -0700
files: src/StokesFACOps.h src/StokesFACOps/computeCompositeResidualOnLevel.C src/StokesFACOps/residual_3D.C src/StokesFACOps/smoothError.C src/StokesFACOps/smooth_Tackley_3D.C src/StokesFACOps/smooth_V_3D.C src/dRm_dv.h
description:
Make it work in 3D
diff -r a01e442d3f2c -r 8e656c4c1517 src/StokesFACOps.h
--- a/src/StokesFACOps.h Tue May 01 15:13:55 2012 -0700
+++ b/src/StokesFACOps.h Tue May 01 15:15:15 2012 -0700
@@ -493,6 +493,13 @@ private:
int num_sweeps,
double residual_tolerance = -1.0);
+ void smooth_Tackley_3D
+ (SAMRAIVectorReal<double>& solution,
+ const SAMRAIVectorReal<double>& residual,
+ int ln,
+ int num_sweeps,
+ double residual_tolerance = -1.0);
+
void smooth_V_2D(const int &axis,
const hier::Box &pbox,
tbox::Pointer<geom::CartesianPatchGeometry> &geom,
@@ -507,6 +514,19 @@ void smooth_V_2D(const int &axis,
pdat::CellData<double> &cell_moduli,
pdat::NodeData<double> &edge_moduli,
const double &theta_momentum);
+
+void smooth_V_3D(const int &ix,
+ const hier::Box &pbox,
+ tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ pdat::SideData<double> &v,
+ pdat::SideData<double> &v_rhs,
+ pdat::CellData<double> &cell_viscosity,
+ pdat::EdgeData<double> &edge_viscosity,
+ const pdat::CellIndex ¢er,
+ const double Dx[3],
+ const double &theta_momentum,
+ const hier::Index pp[3],
+ double &maxres);
/* The mixed derivative of the stress. We have to use a template
because 2D uses Node's for the edge moduli, while 3D uses
diff -r a01e442d3f2c -r 8e656c4c1517 src/StokesFACOps/computeCompositeResidualOnLevel.C
--- a/src/StokesFACOps/computeCompositeResidualOnLevel.C Tue May 01 15:13:55 2012 -0700
+++ b/src/StokesFACOps/computeCompositeResidualOnLevel.C Tue May 01 15:15:15 2012 -0700
@@ -93,6 +93,10 @@ void SAMRAI::solv::StokesFACOps::compute
residual_2D(*p_ptr,*v_ptr,*cell_moduli_ptr,*p_rhs_ptr,*v_rhs_ptr,
*p_resid_ptr,*v_resid_ptr,*patch,pbox,*geom);
break;
+ case 3:
+ residual_3D(*p_ptr,*v_ptr,*cell_moduli_ptr,*p_rhs_ptr,*v_rhs_ptr,
+ *p_resid_ptr,*v_resid_ptr,*patch,pbox,*geom);
+ break;
default:
abort();
}
diff -r a01e442d3f2c -r 8e656c4c1517 src/StokesFACOps/residual_3D.C
--- a/src/StokesFACOps/residual_3D.C Tue May 01 15:13:55 2012 -0700
+++ b/src/StokesFACOps/residual_3D.C Tue May 01 15:15:15 2012 -0700
@@ -4,7 +4,7 @@ void SAMRAI::solv::StokesFACOps::residua
void SAMRAI::solv::StokesFACOps::residual_3D
(pdat::CellData<double> &p,
pdat::SideData<double> &v,
- pdat::CellData<double> &cell_viscosity,
+ pdat::CellData<double> &cell_moduli,
pdat::CellData<double> &p_rhs,
pdat::SideData<double> &v_rhs,
pdat::CellData<double> &p_resid,
@@ -14,8 +14,8 @@ void SAMRAI::solv::StokesFACOps::residua
const geom::CartesianPatchGeometry &geom)
{
tbox::Pointer<pdat::EdgeData<double> >
- edge_viscosity_ptr = patch.getPatchData(edge_viscosity_id);
- pdat::EdgeData<double> &edge_viscosity(*edge_viscosity_ptr);
+ edge_moduli_ptr = patch.getPatchData(edge_moduli_id);
+ pdat::EdgeData<double> &edge_moduli(*edge_moduli_ptr);
const double *Dx = geom.getDx();
const hier::Index ip(1,0,0), jp(0,1,0), kp(0,0,1);
@@ -38,15 +38,7 @@ void SAMRAI::solv::StokesFACOps::residua
if(center[0]!=pbox.upper(0) && center[1]!=pbox.upper(1)
&& center[2]!=pbox.upper(2))
{
- const pdat::SideIndex
- x(center,0,pdat::SideIndex::Lower),
- y(center,1,pdat::SideIndex::Lower),
- z(center,2,pdat::SideIndex::Lower);
-
- double dvx_dx=(v(x+ip) - v(x))/Dx[0];
- double dvy_dy=(v(y+jp) - v(y))/Dx[1];
- double dvz_dz=(v(z+kp) - v(z))/Dx[2];
- p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy - dvz_dz;
+ p_resid(center)=0;
}
for(int ix=0;ix<3;++ix)
@@ -70,7 +62,7 @@ void SAMRAI::solv::StokesFACOps::residua
else
{
v_resid(x)=v_rhs(x)
- - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
+ - v_operator_3D(v,cell_moduli,edge_moduli,
center,edge_y,edge_z,x,y,z,
pp[ix],pp[iy],pp[iz],Dx[ix],Dx[iy],Dx[iz]);
}
diff -r a01e442d3f2c -r 8e656c4c1517 src/StokesFACOps/smoothError.C
--- a/src/StokesFACOps/smoothError.C Tue May 01 15:13:55 2012 -0700
+++ b/src/StokesFACOps/smoothError.C Tue May 01 15:15:15 2012 -0700
@@ -57,6 +57,9 @@ void SAMRAI::solv::StokesFACOps::smoothE
if(d_dim.getValue()==2)
smooth_Tackley_2D(data,residual,ln,num_sweeps,
d_residual_tolerance_during_smoothing);
+ else if(d_dim.getValue()==3)
+ smooth_Tackley_3D(data,residual,ln,num_sweeps,
+ d_residual_tolerance_during_smoothing);
else
TBOX_ERROR(d_object_name << ": Invalid dimension in StokesFACOps.");
} else {
diff -r a01e442d3f2c -r 8e656c4c1517 src/StokesFACOps/smooth_Tackley_3D.C
--- a/src/StokesFACOps/smooth_Tackley_3D.C Tue May 01 15:13:55 2012 -0700
+++ b/src/StokesFACOps/smooth_Tackley_3D.C Tue May 01 15:15:15 2012 -0700
@@ -49,7 +49,6 @@ void SAMRAI::solv::StokesFACOps::smooth_
}
double theta_momentum=0.7;
- double theta_continuity=1.0;
/*
* Smooth the number of sweeps specified or until
@@ -93,11 +92,11 @@ void SAMRAI::solv::StokesFACOps::smooth_
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);
+ = patch->getPatchData(cell_moduli_id);
+ pdat::CellData<double> &cell_moduli(*cell_visc_ptr);
tbox::Pointer<pdat::EdgeData<double> > edge_visc_ptr
- = patch->getPatchData(edge_viscosity_id);
- pdat::EdgeData<double> &edge_viscosity(*edge_visc_ptr);
+ = patch->getPatchData(edge_moduli_id);
+ pdat::EdgeData<double> &edge_moduli(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -115,135 +114,14 @@ void SAMRAI::solv::StokesFACOps::smooth_
pdat::CellIndex center(hier::Index(i,j,k));
/* Update v */
- smooth_V_3D(ix,pbox,geom,p,v,v_rhs,cell_viscosity,
- edge_viscosity,center,
+ smooth_V_3D(ix,pbox,geom,v,v_rhs,cell_moduli,
+ edge_moduli,center,
Dx,theta_momentum,pp,maxres);
}
}
}
set_boundaries(invalid_id,v_id,level,true);
}
-
- /* p sweep
- No need for red-black, because dp does not depend on
- the pressure. */
- xeqScheduleGhostFillNoCoarse(invalid_id,v_id,ln);
-
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
- {
- tbox::Pointer<hier::Patch> patch = *pi;
-
- 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::EdgeData<double> > edge_visc_ptr
- = patch->getPatchData(edge_viscosity_id);
- pdat::EdgeData<double> &edge_viscosity(*edge_visc_ptr);
-
- hier::Box pbox=patch->getBox();
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = patch->getPatchGeometry();
- const double *Dx = geom->getDx();
-
- for(pdat::CellIterator ci(pbox); ci; ci++)
- {
- pdat::CellIndex center(*ci);
-
- double delta_R_continuity=p_rhs(center);
- for(int ix=0;ix<3;++ix)
- {
- const pdat::SideIndex x(center,ix,pdat::SideIndex::Lower);
- delta_R_continuity-=(v(x+pp[ix]) - v(x))/Dx[ix];;
- }
-
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,std::fabs(delta_R_continuity));
-
- dp(center)=delta_R_continuity*theta_continuity
- /dRc_dp_3D(pbox,center,cell_viscosity,edge_viscosity,v,Dx,pp);
- p(center)+=dp(center);
- }
- }
- set_boundaries(p_id,invalid_id,level,true);
-
- /* fix v sweep */
- xeqScheduleGhostFillNoCoarse(dp_id,invalid_id,ln);
-
- for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
- {
- tbox::Pointer<hier::Patch> patch = *pi;
-
- tbox::Pointer<pdat::CellData<double> > dp_ptr =
- patch->getPatchData(dp_id);
- pdat::CellData<double> &dp(*dp_ptr);
-
- tbox::Pointer<pdat::SideData<double> > v_ptr =
- patch->getPatchData(v_id);
- pdat::SideData<double> &v(*v_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::EdgeData<double> > edge_visc_ptr
- = patch->getPatchData(edge_viscosity_id);
- pdat::EdgeData<double> &edge_viscosity(*edge_visc_ptr);
-
- hier::Box pbox=patch->getBox();
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = patch->getPatchGeometry();
- const double *Dx=geom->getDx();
-
- pbox.growUpper(hier::IntVector::getOne(d_dim));
-
- for(pdat::CellIterator ci(pbox); ci; ci++)
- {
- pdat::CellIndex center(*ci);
-
- /* Update v */
- for(int ix=0;ix<3;++ix)
- {
- const int iy((ix+1)%3), iz((ix+2)%3);
- if(center[iy]<pbox.upper(iy) && center[iz]<pbox.upper(iz))
- {
- const pdat::SideIndex x(center,ix,pdat::SideIndex::Lower);
- const pdat::EdgeIndex
- edge_y(center,iy,pdat::EdgeIndex::LowerLeft),
- edge_z(center,iz,pdat::EdgeIndex::LowerLeft);
-
- if(!((center[ix]==pbox.lower(ix)
- && v(x-pp[ix])==boundary_value)
- || (center[ix]==pbox.upper(ix)
- && v(x+pp[ix])==boundary_value)))
- v(x)+=(dp(center) - dp(center-pp[ix]))
- /(Dx[ix]*dRm_dv_3D(cell_viscosity,edge_viscosity,center,
- center-pp[ix],edge_y+pp[iz],edge_y,
- edge_z+pp[iy],edge_z,
- Dx[ix],Dx[iy],Dx[iz]));
- }
- }
- }
- }
- /* This is probably not necessary, since everyone always makes
- sure that everything is set before use. */
- set_boundaries(invalid_id,v_id,level,true);
// if (residual_tolerance >= 0.0) {
/*
diff -r a01e442d3f2c -r 8e656c4c1517 src/StokesFACOps/smooth_V_3D.C
--- a/src/StokesFACOps/smooth_V_3D.C Tue May 01 15:13:55 2012 -0700
+++ b/src/StokesFACOps/smooth_V_3D.C Tue May 01 15:15:15 2012 -0700
@@ -11,7 +11,6 @@ void SAMRAI::solv::StokesFACOps::smooth_
(const int &ix,
const hier::Box &pbox,
tbox::Pointer<geom::CartesianPatchGeometry> &geom,
- pdat::CellData<double> &p,
pdat::SideData<double> &v,
pdat::SideData<double> &v_rhs,
pdat::CellData<double> &cell_viscosity,
@@ -56,7 +55,7 @@ void SAMRAI::solv::StokesFACOps::smooth_
Dx[ix],Dx[iy],Dx[iz]);
double delta_Rx=v_rhs(x)
- - v_operator_3D(v,p,cell_viscosity,edge_viscosity,center,edge_y,edge_z,
+ - v_operator_3D(v,cell_viscosity,edge_viscosity,center,edge_y,edge_z,
x,y,z,pp[ix],pp[iy],pp[iz],Dx[ix],Dx[iy],Dx[iz]);
/* No scaling here, though there should be. */
diff -r a01e442d3f2c -r 8e656c4c1517 src/dRm_dv.h
--- a/src/dRm_dv.h Tue May 01 15:13:55 2012 -0700
+++ b/src/dRm_dv.h Tue May 01 15:15:15 2012 -0700
@@ -23,4 +23,20 @@ double dRm_dv_2D(SAMRAI::pdat::CellData<
-(edge_moduli(up_e,1) + edge_moduli(center_e,1))/(dy*dy);
}
+inline double dRm_dv_3D(SAMRAI::pdat::CellData<double> &cell_moduli,
+ SAMRAI::pdat::EdgeData<double> &edge_moduli,
+ const SAMRAI::pdat::CellIndex ¢er,
+ const SAMRAI::pdat::CellIndex &left,
+ const SAMRAI::pdat::EdgeIndex &front_y,
+ const SAMRAI::pdat::EdgeIndex ¢er_y,
+ const SAMRAI::pdat::EdgeIndex &up_z,
+ const SAMRAI::pdat::EdgeIndex ¢er_z,
+ const double &dx,
+ const double &dy,
+ const double &dz)
+{
+ return dRm_dv_2D(cell_moduli,edge_moduli,center,left,front_y,center_y,dx,dy)
+ - (edge_moduli(up_z,1) + edge_moduli(center_z,1))/(dz*dz);
+}
+
#endif
More information about the CIG-COMMITS
mailing list