[cig-commits] commit: two elastic moduli, still Stokes
Mercurial
hg at geodynamics.org
Thu Jun 7 13:35:34 PDT 2012
changeset: 255:567a621add59
parent: 246:f089bf615532
user: Sylvain Barbot <sbarbot at caltech.edu>
date: Tue Apr 10 11:15:07 2012 -0700
files: input/sinker.input src/FACStokes.h src/FACStokes/FACStokes.C src/FACStokes/fix_moduli.C src/FACStokes/fix_viscosity.C src/FACStokes/initializeLevelData.C src/FACStokes/packDerivedDataIntoDoubleBuffer.C src/FACStokes/setupPlotter.C src/FACStokes/solveStokes.C src/P_MDPI_Refine.C src/P_MDPI_Refine.h src/Resid_Coarsen.C src/Resid_Coarsen.h src/StokesFACOps.h src/StokesFACOps/StokesFACOps.C src/StokesFACOps/computeCompositeResidualOnLevel.C src/StokesFACOps/residual_2D.C src/StokesFACOps/smoothError.C src/StokesFACOps/smooth_Gerya.C src/StokesFACOps/smooth_Tackley_2D.C src/StokesFACOps/smooth_V_2D.C src/StokesFACSolver.h src/StokesFACSolver/initializeSolverState.C src/StokesFACSolver/solveSystem.C src/dRc_dp.h src/dRm_dv.h src/main.C src/set_boundary.C src/viscosity_coarsen.h wscript
description:
two elastic moduli, still Stokes
diff -r f089bf615532 -r 567a621add59 input/sinker.input
--- a/input/sinker.input Thu May 05 13:22:05 2011 -0700
+++ b/input/sinker.input Tue Apr 10 11:15:07 2012 -0700
@@ -29,22 +29,29 @@ FACStokes {
// The inputs for FACStokes is simply the inputs for the individual
// parts owned by the FACStokes class.
adaption_threshold = 1.0e-3
- min_full_refinement_level = 1
+ min_full_refinement_level = 3
- viscosity_ijk= 11, 11
- viscosity_coord_min= -0.001, -0.001
- viscosity_coord_max= 1.001, 1.001
- viscosity_data= 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ lambda_ijk= 11, 11
+ lambda_coord_min= -0.001, -0.001
+ lambda_coord_max= 1.001, 1.001
+ lambda_data= 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1e2, 1e2, 1, 1, 1, 1, 1,
- 1, 1, 1, 1, 1e2, 1e2, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
+
+ mu_ijk= 3, 3
+ mu_coord_min= -0.001, -0.001
+ mu_coord_max= 1.001, 1.001
+ mu_data= 1, 1, 1,
+ 1, 1, 1,
+ 1, 1, 1
v_rhs_ijk= 11, 11
v_rhs_coord_min= -0.001, -0.001
@@ -91,8 +98,8 @@ FACStokes {
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0.03, 0.03, 0, 0, 0, 0, 0,
- 0, 0, 0, 0, 0.03, 0.03, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 3, -3, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 3, -3, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
diff -r f089bf615532 -r 567a621add59 src/FACStokes.h
--- a/src/FACStokes.h Thu May 05 13:22:05 2011 -0700
+++ b/src/FACStokes.h Tue Apr 10 11:15:07 2012 -0700
@@ -141,7 +141,7 @@ namespace SAMRAI {
#endif
private:
- void fix_viscosity();
+ void fix_moduli();
std::string d_object_name;
const tbox::Dimension d_dim;
@@ -185,11 +185,14 @@ namespace SAMRAI {
double d_adaption_threshold;
int min_full_refinement_level;
public:
- int p_id, cell_viscosity_id, edge_viscosity_id, dp_id, p_exact_id,
+ int p_id, cell_moduli_id, edge_moduli_id, dp_id, p_exact_id,
p_rhs_id, v_id, v_rhs_id;
- tbox::Array<double> viscosity, viscosity_xyz_max, viscosity_xyz_min;
- tbox::Array<int> viscosity_ijk;
+ tbox::Array<double> lambda, lambda_xyz_max, lambda_xyz_min;
+ tbox::Array<int> lambda_ijk;
+
+ tbox::Array<double> mu, mu_xyz_max, mu_xyz_min;
+ tbox::Array<int> mu_ijk;
tbox::Array<double> v_rhs, v_rhs_xyz_max, v_rhs_xyz_min;
tbox::Array<int> v_rhs_ijk;
diff -r f089bf615532 -r 567a621add59 src/FACStokes/FACStokes.C
--- a/src/FACStokes/FACStokes.C Thu May 05 13:22:05 2011 -0700
+++ b/src/FACStokes/FACStokes.C Tue Apr 10 11:15:07 2012 -0700
@@ -96,11 +96,12 @@ namespace SAMRAI {
/* ghost cell width is 1 for
stencil widths */);
+ int depth=2;
tbox::Pointer<pdat::CellVariable<double> >
- cell_viscosity_ptr(new pdat::CellVariable<double>(d_dim,
+ cell_moduli_ptr(new pdat::CellVariable<double>(d_dim,
object_name
- + ":cell_viscosity"));
- cell_viscosity_id = vdb->registerVariableAndContext(cell_viscosity_ptr,
+ + ":cell_moduli",depth));
+ cell_moduli_id = vdb->registerVariableAndContext(cell_moduli_ptr,
d_context,
hier::IntVector(d_dim, 1)
/* ghost cell width is
@@ -109,11 +110,11 @@ namespace SAMRAI {
if(dim==2)
{
tbox::Pointer<pdat::NodeVariable<double> >
- edge_viscosity_ptr(new pdat::NodeVariable<double>(d_dim,
+ edge_moduli_ptr(new pdat::NodeVariable<double>(d_dim,
object_name
- + ":edge_viscosity"));
- edge_viscosity_id =
- vdb->registerVariableAndContext(edge_viscosity_ptr,d_context,
+ + ":edge_moduli",depth));
+ edge_moduli_id =
+ vdb->registerVariableAndContext(edge_moduli_ptr,d_context,
hier::IntVector(d_dim,1)
/* ghost cell width is 1 in
case needed */);
@@ -121,11 +122,11 @@ namespace SAMRAI {
else if(dim==3)
{
tbox::Pointer<pdat::EdgeVariable<double> >
- edge_viscosity_ptr(new pdat::EdgeVariable<double>(d_dim,
+ edge_moduli_ptr(new pdat::EdgeVariable<double>(d_dim,
object_name
- + ":edge_viscosity"));
- edge_viscosity_id =
- vdb->registerVariableAndContext(edge_viscosity_ptr,d_context,
+ + ":edge_moduli",depth));
+ edge_moduli_id =
+ vdb->registerVariableAndContext(edge_moduli_ptr,d_context,
hier::IntVector(d_dim,1)
/* ghost cell width is 1 in
case needed */);
@@ -172,14 +173,24 @@ namespace SAMRAI {
min_full_refinement_level
=database->getIntegerWithDefault("min_full_refinement_level",0);
- if(database->keyExists("viscosity_data"))
+ if(database->keyExists("lambda_data"))
{
- viscosity_ijk=database->getIntegerArray("viscosity_ijk");
- viscosity_xyz_min=database->getDoubleArray("viscosity_coord_min");
- viscosity_xyz_max=database->getDoubleArray("viscosity_coord_max");
- viscosity=database->getDoubleArray("viscosity_data");
- check_array_sizes(viscosity_ijk,viscosity_xyz_min,viscosity_xyz_max,
- viscosity,dim,"viscosity");
+ lambda_ijk=database->getIntegerArray("lambda_ijk");
+ lambda_xyz_min=database->getDoubleArray("lambda_coord_min");
+ lambda_xyz_max=database->getDoubleArray("lambda_coord_max");
+ lambda=database->getDoubleArray("lambda_data");
+ check_array_sizes(lambda_ijk,lambda_xyz_min,lambda_xyz_max,
+ lambda,dim,"lambda");
+ }
+
+ if(database->keyExists("mu_data"))
+ {
+ mu_ijk=database->getIntegerArray("mu_ijk");
+ mu_xyz_min=database->getDoubleArray("mu_coord_min");
+ mu_xyz_max=database->getDoubleArray("mu_coord_max");
+ mu=database->getDoubleArray("mu_data");
+ check_array_sizes(mu_ijk,mu_xyz_min,mu_xyz_max,
+ mu,dim,"mu");
}
if(database->keyExists("v_rhs_data"))
diff -r f089bf615532 -r 567a621add59 src/FACStokes/fix_moduli.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/FACStokes/fix_moduli.C Tue Apr 10 11:15:07 2012 -0700
@@ -0,0 +1,138 @@
+#include "FACStokes.h"
+#include "SAMRAI/geom/CartesianGridGeometry.h"
+
+/* Fix the moduli on the coarse grids by coarsening from the finer
+ grids. */
+
+void SAMRAI::FACStokes::fix_moduli()
+{
+ const int ln_max(d_hierarchy->getFinestLevelNumber());
+
+ tbox::Pointer<xfer::CoarsenOperator> cell_moduli_coarsen_operator;
+ tbox::Pointer<xfer::CoarsenAlgorithm> cell_moduli_coarsen_algorithm;
+ tbox::Array<tbox::Pointer<xfer::CoarsenSchedule> >
+ cell_moduli_coarsen_schedules;
+
+ hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
+ tbox::Pointer<geom::CartesianGridGeometry> geometry =
+ d_hierarchy->getGridGeometry();
+ tbox::Pointer<hier::Variable> variable;
+ vdb->mapIndexToVariable(cell_moduli_id, variable);
+ cell_moduli_coarsen_operator =
+ geometry->lookupCoarsenOperator(variable,
+ "CONSERVATIVE_COARSEN");
+ // "CELL_VISCOSITY_COARSEN");
+
+ if (!cell_moduli_coarsen_operator) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot find cell moduli coarsening operator");
+ }
+
+ cell_moduli_coarsen_schedules.resizeArray(ln_max + 1);
+ cell_moduli_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
+ cell_moduli_coarsen_algorithm->
+ registerCoarsen(cell_moduli_id,cell_moduli_id,
+ cell_moduli_coarsen_operator);
+
+ for (int dest_ln = 0; dest_ln < ln_max; ++dest_ln) {
+ cell_moduli_coarsen_schedules[dest_ln] =
+ cell_moduli_coarsen_algorithm->
+ createSchedule(d_hierarchy->getPatchLevel(dest_ln),
+ d_hierarchy->getPatchLevel(dest_ln + 1));
+ if (!cell_moduli_coarsen_schedules[dest_ln]) {
+ TBOX_ERROR(d_object_name
+ << ": Cannot create a coarsen schedule for cell moduli restriction!\n");
+ }
+ }
+
+ for(int dest_ln=ln_max-1; dest_ln>=0; --dest_ln)
+ {
+ xfer::CoarsenAlgorithm coarsener(d_dim);
+ coarsener.registerCoarsen(cell_moduli_id, cell_moduli_id,
+ cell_moduli_coarsen_operator);
+ coarsener.resetSchedule(cell_moduli_coarsen_schedules[dest_ln]);
+ cell_moduli_coarsen_schedules[dest_ln]->coarsenData();
+ cell_moduli_coarsen_algorithm->
+ resetSchedule(cell_moduli_coarsen_schedules[dest_ln]);
+ }
+
+ cell_moduli_coarsen_algorithm.setNull();
+ cell_moduli_coarsen_schedules.setNull();
+
+ /* Compute edge_moduli by averaging the cell moduli. */
+
+ hier::Index ip(hier::Index::getZeroIndex(d_dim)), jp(ip), kp(ip);
+ ip[0]=1;
+ jp[1]=1;
+ if(d_dim.getValue()>2)
+ kp[2]=1;
+ hier::Index pp[]={ip,jp,kp};
+
+ for (int ln = 0; ln <= d_hierarchy->getFinestLevelNumber(); ++ln)
+{
+ tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+ hier::PatchLevel::Iterator i_p(*level);
+ for ( ; i_p; i_p++)
+ {
+ tbox::Pointer<hier::Patch> patch = *i_p;
+ tbox::Pointer<pdat::CellData<double> >
+ cell_moduli_ptr = patch->getPatchData(cell_moduli_id);
+ pdat::CellData<double> &cell_moduli(*cell_moduli_ptr);
+ if(2==d_dim.getValue())
+ {
+ tbox::Pointer<pdat::NodeData<double> >
+ edge_moduli_ptr = patch->getPatchData(edge_moduli_id);
+ pdat::NodeData<double> &edge_moduli(*edge_moduli_ptr);
+
+ for(pdat::NodeIterator ni(edge_moduli.getBox()); ni; ni++)
+ {
+ for (int m=0;m<2;++m)
+ {
+ pdat::NodeIndex e=ni();
+ pdat::CellIndex c(e);
+ cell_moduli(c,m);
+ cell_moduli(c-ip,m);
+ cell_moduli(c-jp,m);
+ cell_moduli(c-ip-jp,m);
+ edge_moduli(e,m)=
+ pow(cell_moduli(c,m)*cell_moduli(c-ip,m)
+ *cell_moduli(c-jp,m)*cell_moduli(c-ip-jp,m),0.25);
+ }
+ }
+ }
+ else
+ {
+ tbox::Pointer<pdat::EdgeData<double> >
+ edge_moduli_ptr = patch->getPatchData(edge_moduli_id);
+ pdat::EdgeData<double> &edge_moduli(*edge_moduli_ptr);
+ for(int axis=0;axis<3;++axis)
+ {
+ const int axis2((axis+1)%3), axis3((axis+2)%3);
+ for(pdat::EdgeIterator ni(edge_moduli.getBox(),axis); ni; ni++)
+ {
+ pdat::EdgeIndex e=ni();
+ pdat::CellIndex c(e);
+ for (int m=0;m<2;++m)
+ {
+ edge_moduli(e,m)=
+ pow(cell_moduli(c,m)*cell_moduli(c-pp[axis2],m)
+ *cell_moduli(c-pp[axis3],m)
+ *cell_moduli(c-pp[axis2]-pp[axis3],m),0.25);
+ }
+ }
+ }
+ }
+ }
+
+ /* Ghost fill */
+ xfer::RefineAlgorithm refiner(d_dim);
+ refiner.registerRefine(edge_moduli_id,edge_moduli_id,
+ edge_moduli_id,
+ tbox::Pointer<xfer::RefineOperator>(0));
+
+ tbox::Pointer<xfer::RefineSchedule> schedule=
+ refiner.createSchedule(d_hierarchy->getPatchLevel(ln));
+
+ schedule->fillData(0.0,false);
+ }
+}
diff -r f089bf615532 -r 567a621add59 src/FACStokes/fix_viscosity.C
--- a/src/FACStokes/fix_viscosity.C Thu May 05 13:22:05 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,132 +0,0 @@
-#include "FACStokes.h"
-#include "SAMRAI/geom/CartesianGridGeometry.h"
-
-/* Fix the viscosity on the coarse grids by coarsening from the finer
- grids. */
-
-void SAMRAI::FACStokes::fix_viscosity()
-{
- const int ln_max(d_hierarchy->getFinestLevelNumber());
-
- tbox::Pointer<xfer::CoarsenOperator> cell_viscosity_coarsen_operator;
- tbox::Pointer<xfer::CoarsenAlgorithm> cell_viscosity_coarsen_algorithm;
- tbox::Array<tbox::Pointer<xfer::CoarsenSchedule> >
- cell_viscosity_coarsen_schedules;
-
- hier::VariableDatabase* vdb = hier::VariableDatabase::getDatabase();
- tbox::Pointer<geom::CartesianGridGeometry> geometry =
- d_hierarchy->getGridGeometry();
- tbox::Pointer<hier::Variable> variable;
- vdb->mapIndexToVariable(cell_viscosity_id, variable);
- cell_viscosity_coarsen_operator =
- geometry->lookupCoarsenOperator(variable,
- "CONSERVATIVE_COARSEN");
- // "CELL_VISCOSITY_COARSEN");
-
- if (!cell_viscosity_coarsen_operator) {
- TBOX_ERROR(d_object_name
- << ": Cannot find cell viscosity coarsening operator");
- }
-
- cell_viscosity_coarsen_schedules.resizeArray(ln_max + 1);
- cell_viscosity_coarsen_algorithm = new xfer::CoarsenAlgorithm(d_dim);
- cell_viscosity_coarsen_algorithm->
- registerCoarsen(cell_viscosity_id,cell_viscosity_id,
- cell_viscosity_coarsen_operator);
-
- for (int dest_ln = 0; dest_ln < ln_max; ++dest_ln) {
- cell_viscosity_coarsen_schedules[dest_ln] =
- cell_viscosity_coarsen_algorithm->
- createSchedule(d_hierarchy->getPatchLevel(dest_ln),
- d_hierarchy->getPatchLevel(dest_ln + 1));
- if (!cell_viscosity_coarsen_schedules[dest_ln]) {
- TBOX_ERROR(d_object_name
- << ": Cannot create a coarsen schedule for cell viscosity restriction!\n");
- }
- }
-
- for(int dest_ln=ln_max-1; dest_ln>=0; --dest_ln)
- {
- xfer::CoarsenAlgorithm coarsener(d_dim);
- coarsener.registerCoarsen(cell_viscosity_id, cell_viscosity_id,
- cell_viscosity_coarsen_operator);
- coarsener.resetSchedule(cell_viscosity_coarsen_schedules[dest_ln]);
- cell_viscosity_coarsen_schedules[dest_ln]->coarsenData();
- cell_viscosity_coarsen_algorithm->
- resetSchedule(cell_viscosity_coarsen_schedules[dest_ln]);
- }
-
- cell_viscosity_coarsen_algorithm.setNull();
- cell_viscosity_coarsen_schedules.setNull();
-
- /* Compute edge_viscosity by averaging the cell viscosities. */
-
- hier::Index ip(hier::Index::getZeroIndex(d_dim)), jp(ip), kp(ip);
- ip[0]=1;
- jp[1]=1;
- if(d_dim.getValue()>2)
- kp[2]=1;
- hier::Index pp[]={ip,jp,kp};
-
- for (int ln = 0; ln <= d_hierarchy->getFinestLevelNumber(); ++ln)
- {
- tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
- hier::PatchLevel::Iterator i_p(*level);
- for ( ; i_p; i_p++)
- {
- tbox::Pointer<hier::Patch> patch = *i_p;
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity_ptr = patch->getPatchData(cell_viscosity_id);
- pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
- if(d_dim.getValue()==2)
- {
- tbox::Pointer<pdat::NodeData<double> >
- edge_viscosity_ptr = patch->getPatchData(edge_viscosity_id);
- pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
-
- for(pdat::NodeIterator ni(edge_viscosity.getBox()); ni; ni++)
- {
- pdat::NodeIndex e=ni();
- pdat::CellIndex c(e);
- cell_viscosity(c);
- cell_viscosity(c-ip);
- cell_viscosity(c-jp);
- cell_viscosity(c-ip-jp);
- edge_viscosity(e)=
- pow(cell_viscosity(c)*cell_viscosity(c-ip)
- *cell_viscosity(c-jp)*cell_viscosity(c-ip-jp),0.25);
- }
- }
- else
- {
- tbox::Pointer<pdat::EdgeData<double> >
- edge_viscosity_ptr = patch->getPatchData(edge_viscosity_id);
- pdat::EdgeData<double> &edge_viscosity(*edge_viscosity_ptr);
- for(int axis=0;axis<3;++axis)
- {
- const int axis2((axis+1)%3), axis3((axis+2)%3);
- for(pdat::EdgeIterator ni(edge_viscosity.getBox(),axis); ni; ni++)
- {
- pdat::EdgeIndex e=ni();
- pdat::CellIndex c(e);
- edge_viscosity(e)=
- pow(cell_viscosity(c)*cell_viscosity(c-pp[axis2])
- *cell_viscosity(c-pp[axis3])
- *cell_viscosity(c-pp[axis2]-pp[axis3]),0.25);
- }
- }
- }
- }
-
- /* Ghost fill */
- xfer::RefineAlgorithm refiner(d_dim);
- refiner.registerRefine(edge_viscosity_id,edge_viscosity_id,
- edge_viscosity_id,
- tbox::Pointer<xfer::RefineOperator>(0));
-
- tbox::Pointer<xfer::RefineSchedule> schedule=
- refiner.createSchedule(d_hierarchy->getPatchLevel(ln));
-
- schedule->fillData(0.0,false);
- }
-}
diff -r f089bf615532 -r 567a621add59 src/FACStokes/initializeLevelData.C
--- a/src/FACStokes/initializeLevelData.C Thu May 05 13:22:05 2011 -0700
+++ b/src/FACStokes/initializeLevelData.C Tue Apr 10 11:15:07 2012 -0700
@@ -37,8 +37,8 @@ void SAMRAI::FACStokes::initializeLevelD
if (allocate_data) {
level->allocatePatchData(p_id);
- level->allocatePatchData(cell_viscosity_id);
- level->allocatePatchData(edge_viscosity_id);
+ level->allocatePatchData(cell_moduli_id);
+ level->allocatePatchData(edge_moduli_id);
level->allocatePatchData(dp_id);
level->allocatePatchData(p_rhs_id);
level->allocatePatchData(p_exact_id);
@@ -61,29 +61,50 @@ void SAMRAI::FACStokes::initializeLevelD
geom = patch->getPatchGeometry();
const double *dx=geom->getDx();
- /* Initialize cell viscosity */
- tbox::Pointer<pdat::CellData<double> > cell_viscosity =
- patch->getPatchData(cell_viscosity_id);
+ /* Initialize cell moduli */
+ tbox::Pointer<pdat::CellData<double> > cell_moduli =
+ patch->getPatchData(cell_moduli_id);
- hier::Box cell_visc_box = cell_viscosity->getBox();
- for(pdat::CellIterator ci(cell_viscosity->getGhostBox()); ci; ci++)
+ hier::Box cell_moduli_box = cell_moduli->getBox();
+
+ for(pdat::CellIterator ci(cell_moduli->getGhostBox()); ci; ci++)
{
pdat::CellIndex c=ci();
double xyz[dim];
for(int d=0;d<dim;++d)
xyz[d]=geom->getXLower()[d]
- + dx[d]*(c[d]-cell_visc_box.lower()[d] + 0.5);
+ + dx[d]*(c[d]-cell_moduli_box.lower()[d] + 0.5);
int ijk(0), factor(1);
for(int d=0;d<dim;++d)
{
- int i=static_cast<int>(xyz[d]*(viscosity_ijk[d]-1)
- /(viscosity_xyz_max[d]-viscosity_xyz_min[d]));
- i=std::max(0,std::min(viscosity_ijk[d]-1,i));
+ int i=static_cast<int>(xyz[d]*(lambda_ijk[d]-1)
+ /(lambda_xyz_max[d]-lambda_xyz_min[d]));
+ i=std::max(0,std::min(lambda_ijk[d]-1,i));
ijk+=i*factor;
- factor*=viscosity_ijk[d];
+ factor*=lambda_ijk[d];
}
- (*cell_viscosity)(c)=viscosity[ijk];
+ (*cell_moduli)(c,0)=lambda[ijk];
+ }
+
+ for(pdat::CellIterator ci(cell_moduli->getGhostBox()); ci; ci++)
+ {
+ pdat::CellIndex c=ci();
+ double xyz[dim];
+ for(int d=0;d<dim;++d)
+ xyz[d]=geom->getXLower()[d]
+ + dx[d]*(c[d]-cell_moduli_box.lower()[d] + 0.5);
+
+ int ijk(0), factor(1);
+ for(int d=0;d<dim;++d)
+ {
+ int i=static_cast<int>(xyz[d]*(mu_ijk[d]-1)
+ /(mu_xyz_max[d]-mu_xyz_min[d]));
+ i=std::max(0,std::min(mu_ijk[d]-1,i));
+ ijk+=i*factor;
+ factor*=mu_ijk[d];
+ }
+ (*cell_moduli)(c,1)=mu[ijk];
}
/* I do not think this is actually necessary. */
diff -r f089bf615532 -r 567a621add59 src/FACStokes/packDerivedDataIntoDoubleBuffer.C
--- a/src/FACStokes/packDerivedDataIntoDoubleBuffer.C Thu May 05 13:22:05 2011 -0700
+++ b/src/FACStokes/packDerivedDataIntoDoubleBuffer.C Tue Apr 10 11:15:07 2012 -0700
@@ -35,30 +35,48 @@ namespace SAMRAI {
variable_name,
int depth_id) const
{
- (void)depth_id;
+ pdat::CellData<double>::Iterator icell(region);
+ const hier::Index ip(1,0), jp(0,1);
- pdat::CellData<double>::Iterator icell(region);
+ tbox::Pointer<pdat::SideData<double> > v_ptr;
+ if (variable_name == "Displacement") {
+ v_ptr = patch.getPatchData(v_id);
+ }
+ else if ("Equivalent body force" == variable_name)
+ {
+ v_ptr = patch.getPatchData(v_rhs_id);
+ }
+ else
+ {
+ // Did not register this name.
+ TBOX_ERROR(
+ "Unregistered variable name '" << variable_name << "' in\n"
+ <<
+ "FACStokesX::packDerivedDataIntoDoubleBuffer");
+ }
- if (variable_name == "Error") {
- tbox::Pointer<pdat::CellData<double> > current_solution_ =
- patch.getPatchData(p_id);
- tbox::Pointer<pdat::CellData<double> > exact_solution_ =
- patch.getPatchData(p_exact_id);
- pdat::CellData<double>& current_solution = *current_solution_;
- pdat::CellData<double>& exact_solution = *exact_solution_;
- for ( ; icell; icell++) {
- double diff = (current_solution(*icell) - exact_solution(*icell));
- *buffer = diff;
- buffer = buffer + 1;
- }
- } else {
- // Did not register this name.
- TBOX_ERROR(
- "Unregistered variable name '" << variable_name << "' in\n"
- <<
- "FACStokesX::packDerivedDataIntoDoubleBuffer");
+ pdat::SideData<double>& v = *v_ptr;
+ for ( ; icell; icell++) {
- }
+ pdat::CellIndex center(*icell);
+ const pdat::SideIndex
+ x(center,0,pdat::SideIndex::Lower),
+ y(center,1,pdat::SideIndex::Lower);
+
+ /* Update p */
+ double vx=(v(x+ip) + v(x))/2.;
+ double vy=(v(y+jp) + v(y))/2.;
+
+ if (0==depth_id)
+ {
+ *buffer = vx;
+ }
+ else
+ {
+ *buffer = vy;
+ }
+ buffer = buffer + 1;
+ }
// Return true if this patch has derived data on it.
// False otherwise.
return true;
diff -r f089bf615532 -r 567a621add59 src/FACStokes/setupPlotter.C
--- a/src/FACStokes/setupPlotter.C Thu May 05 13:22:05 2011 -0700
+++ b/src/FACStokes/setupPlotter.C Tue Apr 10 11:15:07 2012 -0700
@@ -39,15 +39,22 @@ namespace SAMRAI {
plotter.registerPlotQuantity("Pressure",
"SCALAR",
p_id);
- plotter.registerDerivedPlotQuantity("Error",
- "SCALAR",
+ plotter.registerDerivedPlotQuantity("Displacement",
+ "VECTOR",
+ (appu::VisDerivedDataStrategy *)this);
+ plotter.registerDerivedPlotQuantity("Equivalent body force",
+ "VECTOR",
(appu::VisDerivedDataStrategy *)this);
plotter.registerPlotQuantity("Exact solution",
"SCALAR",
p_exact_id);
- plotter.registerPlotQuantity("Cell Viscosity",
+ plotter.registerPlotQuantity("Cell lambda",
"SCALAR",
- cell_viscosity_id);
+ cell_moduli_id,0);
+ // this, below, doesn't seem to work.
+ plotter.registerPlotQuantity("Cell mu",
+ "SCALAR",
+ cell_moduli_id,1);
plotter.registerPlotQuantity("Stokes source",
"SCALAR",
p_rhs_id);
diff -r f089bf615532 -r 567a621add59 src/FACStokes/solveStokes.C
--- a/src/FACStokes/solveStokes.C Thu May 05 13:22:05 2011 -0700
+++ b/src/FACStokes/solveStokes.C Tue Apr 10 11:15:07 2012 -0700
@@ -130,10 +130,10 @@ int SAMRAI::FACStokes::solveStokes()
d_stokes_fac_solver.set_boundaries(p_id,v_id,level,false);
}
- fix_viscosity();
+ fix_moduli();
d_stokes_fac_solver.initializeSolverState
- (p_id,cell_viscosity_id,edge_viscosity_id,dp_id,p_rhs_id,v_id,v_rhs_id,
+ (p_id,cell_moduli_id,edge_moduli_id,dp_id,p_rhs_id,v_id,v_rhs_id,
d_hierarchy,0,d_hierarchy->getFinestLevelNumber());
tbox::plog << "solving..." << std::endl;
diff -r f089bf615532 -r 567a621add59 src/P_MDPI_Refine.C
--- a/src/P_MDPI_Refine.C Thu May 05 13:22:05 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,154 +0,0 @@
-/*************************************************************************
- *
- * 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: Linear refine operator for cell-centered double data on
- * a Cartesian mesh.
- *
- ************************************************************************/
-
-#ifndef included_geom_P_MDPI_Refine_C
-#define included_geom_P_MDPI_Refine_C
-
-#include "P_MDPI_Refine.h"
-
-#include <float.h>
-#include <math.h>
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/hier/Index.h"
-#include "SAMRAI/pdat/CellData.h"
-#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/pdat/NodeData.h"
-#include "SAMRAI/tbox/Utilities.h"
-#include "dRc_dp.h"
-
-void SAMRAI::geom::P_MDPI_Refine::refine(
- hier::Patch& fine,
- const hier::Patch& coarse,
- const int dst_component,
- const int src_component,
- const hier::BoxOverlap& fine_overlap,
- const hier::IntVector& ratio) const
-{
- const pdat::CellOverlap* t_overlap =
- dynamic_cast<const pdat::CellOverlap *>(&fine_overlap);
-
- TBOX_ASSERT(t_overlap != NULL);
-
- const hier::BoxList& boxes = t_overlap->getDestinationBoxList();
- for (hier::BoxList::Iterator b(boxes); b; b++) {
- refine(fine,
- coarse,
- dst_component,
- src_component,
- b(),
- ratio);
- }
-}
-
-void SAMRAI::geom::P_MDPI_Refine::refine(
- hier::Patch& fine_patch,
- const hier::Patch& coarse_patch,
- const int dst_component,
- const int src_component,
- const hier::Box& fine_box,
- const hier::IntVector& ratio) const
-{
- const tbox::Dimension& dimension(getDim());
- const int dim(dimension.getValue());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dimension, fine_patch, coarse_patch,
- fine_box, ratio);
-
- tbox::Pointer<pdat::CellData<double> >
- p_ptr = coarse_patch.getPatchData(src_component);
- pdat::CellData<double> &p(*p_ptr);
- tbox::Pointer<pdat::CellData<double> >
- p_fine_ptr = fine_patch.getPatchData(dst_component);
- pdat::CellData<double> &p_fine(*p_fine_ptr);
- tbox::Pointer<pdat::SideData<double> >
- v_ptr = fine_patch.getPatchData(v_id);
- pdat::SideData<double> &v(*v_ptr);
- tbox::Pointer<pdat::CellData<double> >
- cell_viscosity_ptr = fine_patch.getPatchData(cell_viscosity_id);
- pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
- tbox::Pointer<pdat::NodeData<double> > edge_viscosity2D_ptr;
- tbox::Pointer<pdat::EdgeData<double> > edge_viscosity3D_ptr;
- if(dim==2)
- edge_viscosity2D_ptr = fine_patch.getPatchData(edge_viscosity_id);
- else
- edge_viscosity3D_ptr = fine_patch.getPatchData(edge_viscosity_id);
-
-#ifdef DEBUG_CHECK_ASSERTIONS
- TBOX_ASSERT(!p_ptr.isNull());
- TBOX_ASSERT(!p_fine_ptr.isNull());
- TBOX_ASSERT(p.getDepth() == p_fine.getDepth());
-#endif
-
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = fine_patch.getPatchGeometry();
-
- const double *Dx=geom->getDx();
-
- hier::Box interior(p_fine.getBox());
-
- const hier::Box coarse_box
- (hier::Box::coarsen(fine_box,hier::IntVector::getOne(dimension)*2));
-
- hier::Box cell_box(hier::Index::getZeroIndex(dimension),
- hier::Index::getOneIndex(dimension));
-
- for(pdat::CellIterator ci(coarse_box); ci; ci++)
- {
- pdat::CellIndex coarse(*ci);
- pdat::CellIndex fine(coarse*2);
-
- if(interior.contains(fine))
- {
- pdat::CellData<double>
- dRc_dp(cell_box,1,hier::Index::getZeroIndex(dimension));
- double dRc_dp_total(0);
-
- for(pdat::CellIterator ii(cell_box); ii; ii++)
- {
- pdat::CellIndex c_fine(fine+*ii);
-
- if(dim==2)
- {
- pdat::SideIndex x(c_fine,0,pdat::SideIndex::Lower),
- y(c_fine,1,pdat::SideIndex::Lower);
- dRc_dp(*ii)=dRc_dp_2D(fine_box,c_fine,x,y,cell_viscosity,
- *edge_viscosity2D_ptr,v,Dx[0],Dx[1]);
- }
- else
- {
- const hier::Index ip(1,0,0), jp(0,1,0), kp(0,0,1);
- const hier::Index pp[]={ip,jp,kp};
- dRc_dp(*ii)=dRc_dp_3D(fine_box,c_fine,cell_viscosity,
- *edge_viscosity3D_ptr,v,Dx,pp);
- }
- dRc_dp_total+=dRc_dp(*ii,0);
- }
-
- for(pdat::CellIterator ii(cell_box); ii; ii++)
- {
- pdat::CellIndex c_fine(fine+*ii);
- p_fine(c_fine)=p(coarse)*dRc_dp_total
- /(4*(dim-1)*dRc_dp(*ii));
- }
- }
- else
- {
- /* This should never be used as a real value, so we put in
- * a bad value so that we will notice when it happens */
- for(pdat::CellIterator ii(cell_box); ii; ii++)
- {
- pdat::CellIndex c_fine(fine+*ii);
- if(fine_box.contains(c_fine))
- p_fine(c_fine)=boundary_value;
- }
- }
- }
-}
-#endif
diff -r f089bf615532 -r 567a621add59 src/P_MDPI_Refine.h
--- a/src/P_MDPI_Refine.h Thu May 05 13:22:05 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,144 +0,0 @@
-/*************************************************************************
- *
- * 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: Linear refine operator for cell-centered double data on
- * a Cartesian mesh.
- *
- ************************************************************************/
-
-#ifndef included_geom_P_MDPI_Refine
-#define included_geom_P_MDPI_Refine
-
-#include "SAMRAI/SAMRAI_config.h"
-
-#include "SAMRAI/xfer/RefineOperator.h"
-#include "SAMRAI/hier/Box.h"
-#include "SAMRAI/hier/IntVector.h"
-#include "SAMRAI/hier/Patch.h"
-#include "SAMRAI/tbox/Pointer.h"
-#include "SAMRAI/pdat/CellVariable.h"
-
-#include <string>
-
-namespace SAMRAI {
-namespace geom {
-
-/**
- * Class P_MDPI_Refine implements matrix dependent
- * interpolation for cell-centered double patch data defined over a Cartesian
- * mesh. It should work better for variable viscosity.
- *
- * The findRefineOperator() operator function returns true if the input
- * variable is cell-centered double, and the std::string is "P_MDPI_REFINE".
- *
- * @see xfer::RefineOperator
- */
-
-class P_MDPI_Refine:
- public xfer::RefineOperator
-{
-public:
- /**
- * Uninteresting default constructor.
- */
- explicit P_MDPI_Refine(const tbox::Dimension& dim, const int &v,
- const int &cell_viscosity,
- const int &edge_viscosity):
- xfer::RefineOperator(dim, "P_MDPI_REFINE"),
- v_id(v),
- cell_viscosity_id(cell_viscosity),
- edge_viscosity_id(edge_viscosity)
- {
- d_name_id = "P_MDPI_REFINE";
- }
-
-
- /**
- * Uninteresting virtual destructor.
- */
- virtual ~P_MDPI_Refine(){}
-
- /**
- * Return true if the variable and name std::string match cell-centered
- * double linear interpolation; otherwise, return false.
- */
- bool findRefineOperator(const tbox::Pointer<hier::Variable>& var,
- const std::string& op_name) const
- {
- TBOX_DIM_ASSERT_CHECK_ARGS2(*this, *var);
-
- const tbox::Pointer<pdat::CellVariable<double> > cast_var(var);
- if (!cast_var.isNull() && (op_name == d_name_id)) {
- return true;
- } else {
- return false;
- }
- }
- /**
- * Return name std::string identifier of this refinement operator.
- */
- const std::string& getOperatorName() const
- {
- return d_name_id;
- }
-
- /**
- * The priority of cell-centered double linear interpolation is 0.
- * It will be performed before any user-defined interpolation operations.
- */
- int getOperatorPriority() const
- {
- return 0;
- }
-
- /**
- * The stencil width of the linear interpolation operator is the vector
- * of ones. That is, its stencil extends one cell outside the fine box.
- */
- hier::IntVector getStencilWidth() const
- {
- return hier::IntVector::getOne(getDim());
- }
-
- /**
- * Refine the source component on the coarse patch to the destination
- * component on the fine patch using the cell-centered double linear
- * interpolation operator. Interpolation is performed on the intersection
- * of the destination patch and the boxes contained in fine_overlap
- * It is assumed that the coarse patch contains sufficient data for the
- * stencil width of the refinement operator.
- */
- void refine(hier::Patch& fine,
- const hier::Patch& coarse,
- const int dst_component,
- const int src_component,
- const hier::BoxOverlap& fine_overlap,
- const hier::IntVector& ratio) const;
-
- /**
- * Refine the source component on the coarse patch to the destination
- * component on the fine patch using the cell-centered double linear
- * interpolation operator. Interpolation is performed on the intersection
- * of the destination patch and the fine box. It is assumed that the
- * coarse patch contains sufficient data for the stencil width of the
- * refinement operator. This differs from the above refine() method
- * only in that it operates on a single fine box instead of a BoxOverlap.
- */
- void refine(hier::Patch& fine,
- const hier::Patch& coarse,
- const int dst_component,
- const int src_component,
- const hier::Box& fine_box,
- const hier::IntVector& ratio) const;
-
-private:
- std::string d_name_id;
- int v_id, cell_viscosity_id, edge_viscosity_id;
-};
-
-}
-}
-#endif
diff -r f089bf615532 -r 567a621add59 src/Resid_Coarsen.C
--- a/src/Resid_Coarsen.C Thu May 05 13:22:05 2011 -0700
+++ b/src/Resid_Coarsen.C Tue Apr 10 11:15:07 2012 -0700
@@ -1,11 +1,11 @@
#include "Resid_Coarsen.h"
/**
- * Coarsens using the viscosities as weights. So in 2D
- resid_coarse = (resid(i,j)*viscosity(i,j)
- + resid(i,j+1)*viscosity(i,j+1)
- + resid(i+1,j)*viscosity(i+1,j)
- + resid(i+1,j+1)*viscosity(i+1,j+1))/(4*viscosity_coarse)
+ * Coarsens using the moduli as weights. So in 2D
+ resid_coarse = (resid(i,j)*moduli(i,j)
+ + resid(i,j+1)*moduli(i,j+1)
+ + resid(i+1,j)*moduli(i+1,j)
+ + resid(i+1,j+1)*moduli(i+1,j+1))/(4*moduli_coarse)
*/
void SAMRAI::geom::Resid_Coarsen::coarsen(hier::Patch& coarse,
@@ -25,8 +25,8 @@ void SAMRAI::geom::Resid_Coarsen::coarse
r_ptr = coarse.getPatchData(dst_component);
pdat::CellData<double> &r(*r_ptr);
tbox::Pointer<pdat::CellData<double> >
- cell_viscosity_fine_ptr = fine.getPatchData(cell_viscosity_id);
- pdat::CellData<double> &cell_viscosity_fine(*cell_viscosity_fine_ptr);
+ cell_moduli_fine_ptr = fine.getPatchData(cell_moduli_id);
+ pdat::CellData<double> &cell_moduli_fine(*cell_moduli_fine_ptr);
TBOX_ASSERT(!r_ptr.isNull());
TBOX_ASSERT(!r_fine_ptr.isNull());
@@ -40,14 +40,14 @@ void SAMRAI::geom::Resid_Coarsen::coarse
{
pdat::CellIndex coarse(*ci);
pdat::CellIndex fine(coarse*2);
- double temp(0), viscosity_sum(0);
+ double temp(0), moduli_sum(0);
for(pdat::CellIterator ii(cell_box); ii; ii++)
{
pdat::CellIndex i(*ii);
- temp+=r_fine(fine+i)*cell_viscosity_fine(fine+i);
- viscosity_sum+=cell_viscosity_fine(fine+i);
+ temp+=r_fine(fine+i)*(cell_moduli_fine(fine+i,0)+cell_moduli_fine(fine+i,1));
+ moduli_sum+=(cell_moduli_fine(fine+i,0)+cell_moduli_fine(fine+i,1));
}
- r(coarse)=temp/viscosity_sum;
+ r(coarse)=temp/moduli_sum;
}
}
diff -r f089bf615532 -r 567a621add59 src/Resid_Coarsen.h
--- a/src/Resid_Coarsen.h Thu May 05 13:22:05 2011 -0700
+++ b/src/Resid_Coarsen.h Tue Apr 10 11:15:07 2012 -0700
@@ -9,11 +9,11 @@ namespace geom {
namespace geom {
/**
- * Coarsens using the viscosities as weights. So in 2D
- resid_coarse = (resid(i,j)*viscosity(i,j)
- + resid(i,j+1)*viscosity(i,j+1)
- + resid(i+1,j)*viscosity(i+1,j)
- + resid(i+1,j+1)*viscosity(i+1,j+1))/(4*viscosity_coarse)
+ * Coarsens using the moduli as weights. So in 2D
+ resid_coarse = (resid(i,j)*moduli(i,j)
+ + resid(i,j+1)*moduli(i,j+1)
+ + resid(i+1,j)*moduli(i+1,j)
+ + resid(i+1,j+1)*moduli(i+1,j+1))/(4*moduli_coarse)
* @see xfer::CoarsenOperator
*/
@@ -22,9 +22,9 @@ class Resid_Coarsen:
{
public:
explicit Resid_Coarsen(const tbox::Dimension& dim,
- const int &cell_viscosity):
+ const int &cell_moduli):
xfer::CoarsenOperator(dim, "RESID_COARSEN"),
- cell_viscosity_id(cell_viscosity)
+ cell_moduli_id(cell_moduli)
{
d_name_id = "RESID_COARSEN";
}
@@ -68,7 +68,7 @@ public:
private:
std::string d_name_id;
- const int cell_viscosity_id;
+ const int cell_moduli_id;
};
}
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps.h
--- a/src/StokesFACOps.h Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps.h Tue Apr 10 11:15:07 2012 -0700
@@ -290,12 +290,10 @@ 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 &cell_viscosity,
- const int &edge_viscosity, const int &dp)
+ void set_moduli_id(const int &cell_moduli, const int &edge_moduli)
{
- cell_viscosity_id=cell_viscosity;
- edge_viscosity_id=edge_viscosity;
- dp_id=dp;
+ cell_moduli_id=cell_moduli;
+ edge_moduli_id=edge_moduli;
}
//@}
@@ -408,7 +406,7 @@ public:
void residual_2D
(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,
@@ -420,7 +418,7 @@ public:
void 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,
@@ -488,69 +486,36 @@ private:
* @param residual_tolerance the maximum residual considered to be
* converged
*/
- void
- smooth_Tackley_2D(
+ void smooth_Tackley_2D(
SAMRAIVectorReal<double>& error,
const SAMRAIVectorReal<double>& residual,
int ln,
int num_sweeps,
double residual_tolerance = -1.0);
- void
- smooth_Tackley_3D(
- SAMRAIVectorReal<double>& error,
- const SAMRAIVectorReal<double>& residual,
- int ln,
- int num_sweeps,
- double residual_tolerance = -1.0);
-
- void
- smooth_Gerya(
- SAMRAIVectorReal<double>& error,
- 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,
- const pdat::CellIndex ¢er,
- const hier::Index &ip,
- const hier::Index &jp,
- 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);
-
- void smooth_V_3D
- (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,
- pdat::EdgeData<double> &edge_viscosity,
- const pdat::CellIndex ¢er,
- const double dx[3],
- const double &theta_momentum,
- const hier::Index pp[3],
- double &maxres);
+void smooth_V_2D(const int &axis,
+ const hier::Box &pbox,
+ tbox::Pointer<geom::CartesianPatchGeometry> &geom,
+ const pdat::CellIndex ¢er,
+ const hier::Index &ip,
+ const hier::Index &jp,
+ 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_moduli,
+ pdat::NodeData<double> &edge_moduli,
+ const double &theta_momentum);
/* The mixed derivative of the stress. We have to use a template
- because 2D uses Node's for the edge viscosity, while 3D uses
+ because 2D uses Node's for the edge moduli, while 3D uses
Edge's. Written as if it is dtau_xy_dy. */
template<class E_data, class E_index>
double dtau_mixed(pdat::SideData<double> &v,
- const E_data &edge_viscosity,
+ const E_data &edge_moduli,
const pdat::SideIndex &x,
const pdat::SideIndex &y,
const E_index &edge,
@@ -560,9 +525,9 @@ private:
const double &dy)
{
return ((v(x+jp)-v(x))/(dy*dy)
- + (v(y+jp)-v(y+jp-ip))/(dx*dy)) * edge_viscosity(edge+jp)
+ + (v(y+jp)-v(y+jp-ip))/(dx*dy)) * edge_moduli(edge+jp)
- ((v(x)-v(x-jp))/(dy*dy)
- + (v(y)-v(y-ip))/(dx*dy)) * edge_viscosity(edge);
+ + (v(y)-v(y-ip))/(dx*dy)) * edge_moduli(edge);
}
/* The action of the velocity operator. It is written from the
@@ -571,8 +536,8 @@ private:
double v_operator_2D(pdat::SideData<double> &v,
pdat::CellData<double> &p,
- pdat::CellData<double> &cell_viscosity,
- pdat::NodeData<double> &edge_viscosity,
+ pdat::CellData<double> &cell_moduli,
+ pdat::NodeData<double> &edge_moduli,
const pdat::CellIndex ¢er,
const pdat::NodeIndex &edge,
const pdat::SideIndex &x,
@@ -582,9 +547,9 @@ private:
const double &dx,
const double &dy)
{
- double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_viscosity(center)
- - (v(x)-v(x-ip))*cell_viscosity(center-ip))/(dx*dx);
- double dtau_xy_dy=dtau_mixed(v,edge_viscosity,x,y,edge,ip,jp,dx,dy);
+ double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_moduli(center)
+ - (v(x)-v(x-ip))*cell_moduli(center-ip))/(dx*dx);
+ double dtau_xy_dy=dtau_mixed(v,edge_moduli,x,y,edge,ip,jp,dx,dy);
double dp_dx=(p(center)-p(center-ip))/dx;
return dtau_xx_dx + dtau_xy_dy - dp_dx;
@@ -592,8 +557,8 @@ private:
double v_operator_3D(pdat::SideData<double> &v,
pdat::CellData<double> &p,
- pdat::CellData<double> &cell_viscosity,
- pdat::EdgeData<double> &edge_viscosity,
+ pdat::CellData<double> &cell_moduli,
+ pdat::EdgeData<double> &edge_moduli,
const pdat::CellIndex ¢er,
const pdat::EdgeIndex &edge_y,
const pdat::EdgeIndex &edge_z,
@@ -607,10 +572,10 @@ private:
const double &dy,
const double &dz)
{
- double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_viscosity(center)
- - (v(x)-v(x-ip))*cell_viscosity(center-ip))/(dx*dx);
- double dtau_xy_dy=dtau_mixed(v,edge_viscosity,x,y,edge_z,ip,jp,dx,dy);
- double dtau_xz_dz=dtau_mixed(v,edge_viscosity,x,z,edge_y,ip,kp,dx,dz);
+ double dtau_xx_dx=2*((v(x+ip)-v(x))*cell_moduli(center)
+ - (v(x)-v(x-ip))*cell_moduli(center-ip))/(dx*dx);
+ double dtau_xy_dy=dtau_mixed(v,edge_moduli,x,y,edge_z,ip,jp,dx,dy);
+ double dtau_xz_dz=dtau_mixed(v,edge_moduli,x,z,edge_y,ip,kp,dx,dz);
double dp_dx=(p(center)-p(center-ip))/dx;
return dtau_xx_dx + dtau_xy_dy + dtau_xz_dz - dp_dx;
@@ -906,11 +871,11 @@ private:
double d_residual_tolerance_during_smoothing;
/*!
- * @brief Id of viscosity and dp.
+ * @brief Id of moduli and dp.
*
- * @see set_viscosity_dp_id.
+ * @see set_moduli_dp_id.
*/
- int cell_viscosity_id, edge_viscosity_id, dp_id;
+ int cell_moduli_id, edge_moduli_id, dp_id;
#ifdef HAVE_HYPRE
/*!
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/StokesFACOps.C
--- a/src/StokesFACOps/StokesFACOps.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps/StokesFACOps.C Tue Apr 10 11:15:07 2012 -0700
@@ -60,8 +60,8 @@ namespace SAMRAI {
d_coarse_solver_tolerance(1.e-8),
d_coarse_solver_max_iterations(10),
d_residual_tolerance_during_smoothing(-1.0),
- cell_viscosity_id(invalid_id),
- edge_viscosity_id(invalid_id),
+ cell_moduli_id(invalid_id),
+ edge_moduli_id(invalid_id),
dp_id(invalid_id),
#ifdef HAVE_HYPRE
d_hypre_solver(dim,
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/computeCompositeResidualOnLevel.C
--- a/src/StokesFACOps/computeCompositeResidualOnLevel.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps/computeCompositeResidualOnLevel.C Tue Apr 10 11:15:07 2012 -0700
@@ -72,7 +72,7 @@ void SAMRAI::solv::StokesFACOps::compute
tbox::Pointer<pdat::SideData<double> >
v_ptr = solution.getComponentPatchData(1, *patch);
tbox::Pointer<pdat::CellData<double> >
- cell_viscosity_ptr = patch->getPatchData(cell_viscosity_id);
+ cell_moduli_ptr = patch->getPatchData(cell_moduli_id);
tbox::Pointer<pdat::CellData<double> >
p_rhs_ptr = rhs.getComponentPatchData(0, *patch);
tbox::Pointer<pdat::SideData<double> >
@@ -90,11 +90,7 @@ void SAMRAI::solv::StokesFACOps::compute
switch(d_dim.getValue())
{
case 2:
- residual_2D(*p_ptr,*v_ptr,*cell_viscosity_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_viscosity_ptr,*p_rhs_ptr,*v_rhs_ptr,
+ 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;
default:
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/residual_2D.C
--- a/src/StokesFACOps/residual_2D.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps/residual_2D.C Tue Apr 10 11:15:07 2012 -0700
@@ -4,7 +4,7 @@ void SAMRAI::solv::StokesFACOps::residua
void SAMRAI::solv::StokesFACOps::residual_2D
(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,
@@ -16,8 +16,8 @@ void SAMRAI::solv::StokesFACOps::residua
const hier::Index ip(1,0), jp(0,1);
tbox::Pointer<pdat::NodeData<double> >
- edge_viscosity_ptr = patch.getPatchData(edge_viscosity_id);
- pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
+ edge_moduli_ptr = patch.getPatchData(edge_moduli_id);
+ pdat::NodeData<double> &edge_moduli(*edge_moduli_ptr);
double dx = geom.getDx()[0];
double dy = geom.getDx()[1];
@@ -42,9 +42,7 @@ void SAMRAI::solv::StokesFACOps::residua
/* p */
if(center[0]!=pbox.upper(0) && center[1]!=pbox.upper(1))
{
- double dvx_dx=(v(x+ip) - v(x))/dx;
- double dvy_dy=(v(y+jp) - v(y))/dy;
- p_resid(center)=p_rhs(center) - dvx_dx - dvy_dy;
+ p_resid(center)=0;
}
/* vx */
@@ -59,7 +57,7 @@ void SAMRAI::solv::StokesFACOps::residua
else
{
v_resid(x)=v_rhs(x)
- - v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator_2D(v,p,cell_moduli,edge_moduli,center,
edge,x,y,ip,jp,dx,dy);
}
}
@@ -76,7 +74,7 @@ void SAMRAI::solv::StokesFACOps::residua
else
{
v_resid(y)=v_rhs(y)
- - v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator_2D(v,p,cell_moduli,edge_moduli,center,
edge,y,x,jp,ip,dy,dx);
}
}
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/smoothError.C
--- a/src/StokesFACOps/smoothError.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps/smoothError.C Tue Apr 10 11:15:07 2012 -0700
@@ -53,16 +53,12 @@ void SAMRAI::solv::StokesFACOps::smoothE
{
t_smooth_error->start();
- if (d_smoothing_choice == "Gerya") {
- smooth_Gerya(data,residual,ln,num_sweeps,
- d_residual_tolerance_during_smoothing);
- } else if (d_smoothing_choice == "Tackley") {
+ if (d_smoothing_choice == "Tackley") {
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 {
TBOX_ERROR(d_object_name << ": Bad smoothing choice '"
<< d_smoothing_choice
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/smooth_Gerya.C
--- a/src/StokesFACOps/smooth_Gerya.C Thu May 05 13:22:05 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,184 +0,0 @@
-#include "StokesFACOps.h"
-#include "Constants.h"
-/* Smooth the error using a red-black Gauss-Seidel smoother inspired
- by Introduction to Numerical Geodynamic Modelling, Taras Gerya,
- 2010
-
- This does not give the same answers in serial and parallel, because
- it smooths both vx and vy at the same time. The stencil for the vx
- uses diagonal elements from vy and vice-versa. So the red and
- black updates are not entirely separate. */
-
-void SAMRAI::solv::StokesFACOps::smooth_Gerya
-(SAMRAIVectorReal<double>& solution,
- const SAMRAIVectorReal<double>& residual,
- int ln,
- int num_sweeps,
- double residual_tolerance)
-{
- const int p_id(solution.getComponentDescriptorIndex(0)),
- p_rhs_id(residual.getComponentDescriptorIndex(0)),
- v_id(solution.getComponentDescriptorIndex(1)),
- v_rhs_id(residual.getComponentDescriptorIndex(1));
-
-#ifdef DEBUG_CHECK_ASSERTIONS
- if (solution.getPatchHierarchy() != d_hierarchy
- || residual.getPatchHierarchy() != d_hierarchy)
- {
- TBOX_ERROR(d_object_name << ": Vector hierarchy does not match\n"
- "internal hierarchy.");
- }
-#endif
- tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
-
- /* Only need to sync the rhs once. This sync is needed because
- calculating a new pressure update requires computing in the ghost
- region so that the update for the velocity inside the box will be
- correct. */
- p_refine_patch_strategy.setTargetDataId(p_id);
- v_refine_patch_strategy.setTargetDataId(v_id);
- set_boundaries(p_id,v_id,level,true);
- xeqScheduleGhostFillNoCoarse(p_rhs_id,v_rhs_id,ln);
-
- if (ln > d_ln_min) {
- /*
- * Perform a one-time transfer of data from coarser level,
- * to fill ghost boundaries that will not change through
- * the smoothing loop.
- */
- xeqScheduleGhostFill(p_id, v_id, ln);
- }
-
- double theta_momentum=1.2;
- double theta_continuity=0.3;
-
- /*
- * Smooth the number of sweeps specified or until
- * the convergence is satisfactory.
- */
- double maxres;
- /*
- * 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.
- */
- const hier::Index ip(1,0), jp(0,1);
- bool converged = false;
- for (int sweep=0; sweep < num_sweeps*(1<<(d_ln_max-ln)) && !converged;
- ++sweep)
- {
- maxres=0;
- for(int rb=0;rb<2;++rb)
- {
- 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_ptr =
- patch->getPatchData(p_id);
- pdat::CellData<double> &p(*p_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>
- geom = patch->getPatchGeometry();
- double dx = geom->getDx()[0];
- double dy = geom->getDx()[1];
-
- 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);
-
- ++up[1];
- --down[1];
- ++right[0];
- --left[0];
-
- const pdat::SideIndex
- x(center,0,pdat::SideIndex::Lower),
- y(center,1,pdat::SideIndex::Lower);
- const pdat::NodeIndex
- edge(center,pdat::NodeIndex::LowerLeft);
-
- /* Update p */
- if(j<pbox.upper(1)+1 && i<pbox.upper(0)+1)
- {
- double dvx_dx=(v(x+ip)-v(x))/dx;
- double dvy_dy=(v(y+jp)-v(y))/dy;
-
- double delta_R_continuity=p_rhs(center)-dvx_dx-dvy_dy;
-
- /* No scaling here, though there should be. */
- maxres=std::max(maxres,std::fabs(delta_R_continuity));
-
- p(center)+=cell_viscosity(center)
- *delta_R_continuity*theta_continuity;
- }
- /* Update v */
- if(j<pbox.upper(1)+1)
- {
- smooth_V_2D(0,pbox,geom,center,ip,jp,
- p,v,v_rhs,maxres,dx,dy,cell_viscosity,
- edge_viscosity,theta_momentum);
- }
- if(i<pbox.upper(0)+1)
- {
- smooth_V_2D(1,pbox,geom,center,jp,ip,
- p,v,v_rhs,maxres,dy,dx,cell_viscosity,
- edge_viscosity,theta_momentum);
- }
- }
- }
- }
- set_boundaries(p_id,v_id,level,true);
- }
- // 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"
- << "Gerya " << ln << " " << sweep << " : " << maxres << "\n";
- // }
- }
-}
-
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/smooth_Tackley_2D.C
--- a/src/StokesFACOps/smooth_Tackley_2D.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps/smooth_Tackley_2D.C Tue Apr 10 11:15:07 2012 -0700
@@ -48,8 +48,7 @@ void SAMRAI::solv::StokesFACOps::smooth_
xeqScheduleGhostFill(p_id, v_id, ln);
}
- double theta_momentum=0.7;
- double theta_continuity=1.0;
+ double theta_momentum=1.0;
/*
* Smooth the number of sweeps specified or until
@@ -90,11 +89,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::NodeData<double> > edge_visc_ptr
- = patch->getPatchData(edge_viscosity_id);
- pdat::NodeData<double> &edge_viscosity(*edge_visc_ptr);
+ = patch->getPatchData(edge_moduli_id);
+ pdat::NodeData<double> &edge_moduli(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -114,8 +113,8 @@ void SAMRAI::solv::StokesFACOps::smooth_
/* Update v */
smooth_V_2D(0,pbox,geom,center,ip,jp,
- p,v,v_rhs,maxres,dx,dy,cell_viscosity,
- edge_viscosity,theta_momentum);
+ p,v,v_rhs,maxres,dx,dy,cell_moduli,
+ edge_moduli,theta_momentum);
}
}
}
@@ -144,11 +143,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::NodeData<double> > edge_visc_ptr
- = patch->getPatchData(edge_viscosity_id);
- pdat::NodeData<double> &edge_viscosity(*edge_visc_ptr);
+ = patch->getPatchData(edge_moduli_id);
+ pdat::NodeData<double> &edge_moduli(*edge_visc_ptr);
hier::Box pbox=patch->getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
@@ -168,137 +167,13 @@ void SAMRAI::solv::StokesFACOps::smooth_
/* Update v */
smooth_V_2D(1,pbox,geom,center,jp,ip,
- p,v,v_rhs,maxres,dy,dx,cell_viscosity,
- edge_viscosity,theta_momentum);
+ p,v,v_rhs,maxres,dy,dx,cell_moduli,
+ edge_moduli,theta_momentum);
}
}
}
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::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>
- geom = patch->getPatchGeometry();
- double dx = geom->getDx()[0];
- double dy = geom->getDx()[1];
-
- for(pdat::CellIterator ci(pbox); ci; ci++)
- {
- pdat::CellIndex center(*ci);
- const pdat::SideIndex
- x(center,0,pdat::SideIndex::Lower),
- y(center,1,pdat::SideIndex::Lower);
-
- /* Update p */
- double dvx_dx=(v(x+ip) - v(x))/dx;
- double dvy_dy=(v(y+jp) - v(y))/dy;
-
- double delta_R_continuity=
- 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)=delta_R_continuity*theta_continuity
- /dRc_dp_2D(pbox,center,x,y,cell_viscosity,edge_viscosity,v,dx,dy);
- 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::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>
- geom = patch->getPatchGeometry();
- double dx = geom->getDx()[0];
- double dy = geom->getDx()[1];
-
- pbox.growUpper(hier::IntVector::getOne(d_dim));
-
- for(pdat::CellIterator ci(pbox); ci; ci++)
- {
- pdat::CellIndex center(*ci);
-
- const pdat::SideIndex x(center,0,pdat::SideIndex::Lower),
- y(center,1,pdat::SideIndex::Lower);
- const pdat::NodeIndex edge(center,pdat::NodeIndex::LowerLeft);
-
- /* Update v */
- if(center[1]<pbox.upper(1))
- {
- if(!((center[0]==pbox.lower(0) && v(x-ip)==boundary_value)
- || (center[0]==pbox.upper(0)
- && v(x+ip)==boundary_value)))
- v(x)+=(dp(center) - dp(center-ip))
- /(dx*dRm_dv_2D(cell_viscosity,edge_viscosity,center,
- center-ip,edge+jp,edge,dx,dy));
- }
- if(center[0]<pbox.upper(0))
- {
- if(!((center[1]==pbox.lower(1) && v(y-jp)==boundary_value)
- || (center[1]==pbox.upper(1)
- && v(y+jp)==boundary_value)))
- v(y)+=(dp(center) - dp(center-jp))
- /(dy*dRm_dv_2D(cell_viscosity,edge_viscosity,center,
- center-jp,edge+ip,edge,dy,dx));
- }
- }
- }
- set_boundaries(invalid_id,v_id,level,true);
// if (residual_tolerance >= 0.0) {
/*
diff -r f089bf615532 -r 567a621add59 src/StokesFACOps/smooth_V_2D.C
--- a/src/StokesFACOps/smooth_V_2D.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACOps/smooth_V_2D.C Tue Apr 10 11:15:07 2012 -0700
@@ -20,8 +20,8 @@ void SAMRAI::solv::StokesFACOps::smooth_
double &maxres,
const double &dx,
const double &dy,
- pdat::CellData<double> &cell_viscosity,
- pdat::NodeData<double> &edge_viscosity,
+ pdat::CellData<double> &cell_moduli,
+ pdat::NodeData<double> &edge_moduli,
const double &theta_momentum)
{
const int off_axis=(axis==0) ? 1 : 0;
@@ -53,11 +53,11 @@ void SAMRAI::solv::StokesFACOps::smooth_
dv_upper=v(x+offset) - v(x);
}
- double C_vx=dRm_dv_2D(cell_viscosity,edge_viscosity,center,center-ip,
+ double C_vx=dRm_dv_2D(cell_moduli,edge_moduli,center,center-ip,
edge+jp,edge,dx,dy);
double delta_Rx=v_rhs(x)
- - v_operator_2D(v,p,cell_viscosity,edge_viscosity,center,
+ - v_operator_2D(v,p,cell_moduli,edge_moduli,center,
edge,x,y,ip,jp,dx,dy);
/* No scaling here, though there should be. */
diff -r f089bf615532 -r 567a621add59 src/StokesFACSolver.h
--- a/src/StokesFACSolver.h Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACSolver.h Tue Apr 10 11:15:07 2012 -0700
@@ -179,8 +179,8 @@ public:
bool
solveSystem(
const int p,
- const int cell_viscosity,
- const int edge_viscosity,
+ const int cell_moduli,
+ const int edge_moduli,
const int dp,
const int p_rhs,
const int v,
@@ -509,8 +509,8 @@ public:
*/
void
initializeSolverState(const int p,
- const int cell_viscosity,
- const int edge_viscosity,
+ const int cell_moduli,
+ const int edge_moduli,
const int dp,
const int p_rhs,
const int v,
diff -r f089bf615532 -r 567a621add59 src/StokesFACSolver/initializeSolverState.C
--- a/src/StokesFACSolver/initializeSolverState.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACSolver/initializeSolverState.C Tue Apr 10 11:15:07 2012 -0700
@@ -30,8 +30,8 @@
void SAMRAI::solv::StokesFACSolver::initializeSolverState
(const int p,
- const int cell_viscosity,
- const int edge_viscosity,
+ const int cell_moduli,
+ const int edge_moduli,
const int dp,
const int p_rhs,
const int v,
@@ -97,7 +97,7 @@ void SAMRAI::solv::StokesFACSolver::init
d_ln_max);
}
- d_fac_ops.set_viscosity_dp_id(cell_viscosity,edge_viscosity,dp);
+ d_fac_ops.set_moduli_id(cell_moduli,edge_moduli);
createVectorWrappers(p, p_rhs, v, v_rhs);
diff -r f089bf615532 -r 567a621add59 src/StokesFACSolver/solveSystem.C
--- a/src/StokesFACSolver/solveSystem.C Thu May 05 13:22:05 2011 -0700
+++ b/src/StokesFACSolver/solveSystem.C Tue Apr 10 11:15:07 2012 -0700
@@ -82,8 +82,8 @@ bool SAMRAI::solv::StokesFACSolver::solv
bool SAMRAI::solv::StokesFACSolver::solveSystem
(const int p,
- const int cell_viscosity,
- const int edge_viscosity,
+ const int cell_moduli,
+ const int edge_moduli,
const int dp,
const int p_rhs,
const int v,
@@ -116,7 +116,7 @@ bool SAMRAI::solv::StokesFACSolver::solv
<< "specified.\n");
}
#endif
- initializeSolverState(p, cell_viscosity, edge_viscosity, dp, p_rhs, v, v_rhs,
+ initializeSolverState(p, cell_moduli, edge_moduli, dp, p_rhs, v, v_rhs,
hierarchy, coarse_ln, fine_ln);
bool solver_rval;
diff -r f089bf615532 -r 567a621add59 src/dRc_dp.h
--- a/src/dRc_dp.h Thu May 05 13:22:05 2011 -0700
+++ b/src/dRc_dp.h Tue Apr 10 11:15:07 2012 -0700
@@ -18,8 +18,8 @@ inline double dRc_dp_2D(const SAMRAI::hi
const SAMRAI::pdat::CellIndex ¢er,
const SAMRAI::pdat::SideIndex &x,
const SAMRAI::pdat::SideIndex &y,
- SAMRAI::pdat::CellData<double> &cell_viscosity,
- SAMRAI::pdat::NodeData<double> &edge_viscosity,
+ SAMRAI::pdat::CellData<double> &cell_moduli,
+ SAMRAI::pdat::NodeData<double> &edge_moduli,
SAMRAI::pdat::SideData<double> &v,
const double &dx,
const double &dy)
@@ -37,60 +37,25 @@ inline double dRc_dp_2D(const SAMRAI::hi
double result(0);
if(!(center[0]==pbox.lower(0) && v(x-ip)==boundary_value))
- result+=dRc_dvx_m * dRm_dp_xm/dRm_dv_2D(cell_viscosity,edge_viscosity,
+ result+=dRc_dvx_m * dRm_dp_xm/dRm_dv_2D(cell_moduli,edge_moduli,
center,center-ip,up_e,center_e,
dx,dy);
if(!(center[0]==pbox.upper(0) && v(x+ip*2)==boundary_value))
- result+=dRc_dvx_p * dRm_dp_xp/dRm_dv_2D(cell_viscosity,edge_viscosity,
+ result+=dRc_dvx_p * dRm_dp_xp/dRm_dv_2D(cell_moduli,edge_moduli,
center+ip,center,up_e+ip,
center_e+ip,dx,dy);
if(!(center[1]==pbox.lower(1) && v(y-jp)==boundary_value))
- result+=dRc_dvy_m * dRm_dp_ym/dRm_dv_2D(cell_viscosity,edge_viscosity,
+ result+=dRc_dvy_m * dRm_dp_ym/dRm_dv_2D(cell_moduli,edge_moduli,
center,center-jp,right_e,center_e,
dy,dx);
if(!(center[1]==pbox.upper(1) && v(y+jp*2)==boundary_value))
- result+=dRc_dvy_p * dRm_dp_yp/dRm_dv_2D(cell_viscosity,edge_viscosity,
+ result+=dRc_dvy_p * dRm_dp_yp/dRm_dv_2D(cell_moduli,edge_moduli,
center+jp,center,right_e+jp,
center_e+jp,dy,dx);
return result;
}
-
-inline double dRc_dp_3D(const SAMRAI::hier::Box &pbox,
- const SAMRAI::pdat::CellIndex ¢er,
- SAMRAI::pdat::CellData<double> &cell_viscosity,
- SAMRAI::pdat::EdgeData<double> &edge_viscosity,
- SAMRAI::pdat::SideData<double> &v,
- const double Dx[3],
- const SAMRAI::hier::Index pp[3])
-{
- double result(0);
- for(int ix=0;ix<3;++ix)
- {
- const int iy((ix+1)%3), iz((ix+2)%3);
- const SAMRAI::pdat::SideIndex x(center,ix,SAMRAI::pdat::SideIndex::Lower);
- const SAMRAI::pdat::EdgeIndex
- center_y(center,iy,SAMRAI::pdat::EdgeIndex::LowerLeft),
- front_y(center_y+pp[iz]),
- center_z(center,iz,SAMRAI::pdat::EdgeIndex::LowerLeft),
- up_z(center_z+pp[iy]);
-
- if(!(center[ix]==pbox.lower(ix) && v(x-pp[ix])==boundary_value))
- result+=-1.0/(dRm_dv_3D(cell_viscosity,edge_viscosity,
- center,center-pp[ix],front_y,center_y,
- up_z,center_z,Dx[ix],Dx[iy],Dx[iz])
- * Dx[ix]*Dx[ix]);
-
- if(!(center[ix]==pbox.upper(ix) && v(x+pp[ix]*2)==boundary_value))
- result+=-1.0/(dRm_dv_3D(cell_viscosity,edge_viscosity,
- center+pp[ix],center,front_y+pp[ix],center_y+pp[ix],
- up_z+pp[ix],center_z+pp[ix],
- Dx[ix],Dx[iy],Dx[iz])
- * Dx[ix]*Dx[ix]);
- }
- return result;
-}
#endif
diff -r f089bf615532 -r 567a621add59 src/dRm_dv.h
--- a/src/dRm_dv.h Thu May 05 13:22:05 2011 -0700
+++ b/src/dRm_dv.h Tue Apr 10 11:15:07 2012 -0700
@@ -9,8 +9,8 @@
different values for center etc. to get vy or vx(!center_x). */
template<class E_data, class E_index>
-double dRm_dv_2D(SAMRAI::pdat::CellData<double> &cell_viscosity,
- E_data &edge_viscosity,
+double dRm_dv_2D(SAMRAI::pdat::CellData<double> &cell_moduli,
+ E_data &edge_moduli,
const SAMRAI::pdat::CellIndex ¢er,
const SAMRAI::pdat::CellIndex &left,
const E_index &up_e,
@@ -18,27 +18,8 @@ double dRm_dv_2D(SAMRAI::pdat::CellData<
const double &dx,
const double &dy)
{
- return -2*(cell_viscosity(center) + cell_viscosity(left))/(dx*dx)
- - (edge_viscosity(up_e) + edge_viscosity(center_e))/(dy*dy);
+ return -2*(cell_moduli(center) + cell_moduli(left))/(dx*dx)
+ - (edge_moduli(up_e) + edge_moduli(center_e))/(dy*dy);
}
-inline double dRm_dv_3D(SAMRAI::pdat::CellData<double> &cell_viscosity,
- SAMRAI::pdat::EdgeData<double> &edge_viscosity,
- 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_viscosity,edge_viscosity,center,left,front_y,center_y,dx,dz)
- - (edge_viscosity(up_z) + edge_viscosity(center_z))/(dy*dy);
-}
-
-
-
-
#endif
diff -r f089bf615532 -r 567a621add59 src/main.C
--- a/src/main.C Thu May 05 13:22:05 2011 -0700
+++ b/src/main.C Tue Apr 10 11:15:07 2012 -0700
@@ -35,7 +35,6 @@ using namespace std;
#include "V_Boundary_Refine.h"
#include "V_Coarsen.h"
#include "Resid_Coarsen.h"
-#include "P_MDPI_Refine.h"
#include "FACStokes.h"
@@ -190,14 +189,9 @@ int main(
input_db->getDatabase("FACStokes") :
tbox::Pointer<tbox::Database>(NULL));
- grid_geometry->addSpatialRefineOperator
- (tbox::Pointer<SAMRAI::xfer::RefineOperator>
- (new SAMRAI::geom::P_MDPI_Refine(dim,fac_stokes.v_id,
- fac_stokes.cell_viscosity_id,
- fac_stokes.edge_viscosity_id)));
grid_geometry->addSpatialCoarsenOperator
(tbox::Pointer<SAMRAI::xfer::CoarsenOperator>
- (new SAMRAI::geom::Resid_Coarsen(dim,fac_stokes.cell_viscosity_id)));
+ (new SAMRAI::geom::Resid_Coarsen(dim,fac_stokes.cell_moduli_id)));
/*
* Create the tag-and-initializer, box-generator and load-balancer
diff -r f089bf615532 -r 567a621add59 src/set_boundary.C
--- a/src/set_boundary.C Thu May 05 13:22:05 2011 -0700
+++ b/src/set_boundary.C Tue Apr 10 11:15:07 2012 -0700
@@ -35,65 +35,6 @@ void set_boundary(const SAMRAI::hier::Pa
bool upper_dirichlet[]={true,true,true};
double upper_boundary[]={0,0,0};
-
- if(p_id!=invalid_id)
- {
- tbox::Pointer<pdat::CellData<double> > p_ptr = patch.getPatchData(p_id);
- pdat::CellData<double> &p(*p_ptr);
-
- hier::Box gbox=p.getGhostBox();
-
- for(pdat::CellIterator ci(gbox); ci; ci++)
- {
- pdat::CellIndex center(*ci);
- hier::Index ip(zero), jp(zero), kp(zero);
- hier::Index pp[]={ip,jp,kp};
-
- /* Check if we are at a boundary. If it is a Neumann
- * boundary and not on a corner, then use Neumann
- * conditions. Otherwise use extrapolation */
-
- bool dirichlet_boundary(false);
-
- for(int ix=0;ix<dim;++ix)
- {
- if(center[ix]<pbox.lower(ix)
- && geom->getTouchesRegularBoundary(ix,0))
- {
- pp[ix][ix]=1;
- if(!lower_dirichlet[ix])
- {
- p(center)=-p(center+pp[ix]);
- if(rhs)
- p(center)+=2*p_lower[ix];
- }
- else
- dirichlet_boundary=true;
- }
- else if(center[ix]>pbox.upper(ix)
- && geom->getTouchesRegularBoundary(ix,1))
- {
- pp[ix][ix]=-1;
- if(!upper_dirichlet[ix])
- {
- p(center)=-p(center+pp[ix]);
- if(rhs)
- p(center)+=2*p_upper[ix];
- }
- else
- dirichlet_boundary=true;
- }
- }
-
- /* Actually apply the BC for Dirichlet velocities. */
- if(dirichlet_boundary)
- {
- p(center)=
- 2*p(center+pp[0]+pp[1]+pp[2])-p(center+(pp[0]+pp[1]+pp[2])*2);
- }
- }
- }
-
if(v_id!=invalid_id)
{
diff -r f089bf615532 -r 567a621add59 src/viscosity_coarsen.h
--- a/src/viscosity_coarsen.h Thu May 05 13:22:05 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-#ifndef GAMR_VISCOSITY_COARSEN_H
-#define GAMR_VISCOSITY_COARSEN_H
-
-/* This uses both the cell-centered and node-centered viscosities to
- coarsen the viscosity. In 2D, if you rotate the grid 45 degrees,
- then it becomes a regular node-based grid. So we use a standard
- full-weighted coarsening. There is no refining needed, since
- this is only for viscosity.
-
- Note that coarsening for the node and for the cell are the same
- except for an offset.
-
- Ee------e-------Ee------e-------Ee
- | | | | |
- | c | c | c | c |
- | | | | |
- e-------Ce------e-------Ce------e
- | | | | |
- | c | c | c | c |
- | | | | |
- Ee------e-------Ee------e-------Ee
- | | | | |
- | c | c | c | c |
- | | | | |
- e-------Ce------e-------Ce------e
- | | | | |
- | c | c | c | c |
- | | | | |
- Ee------e-------Ee------e-------Ee
-
-*/
-
-#include "SAMRAI/pdat/CellVariable.h"
-#include "SAMRAI/pdat/NodeVariable.h"
-#include "SAMRAI/pdat/EdgeVariable.h"
-#include "SAMRAI/hier/Index.h"
-
-inline double viscosity_coarsen_2D(const SAMRAI::pdat::CellData<double> &cell,
- const SAMRAI::pdat::NodeData<double> &edge,
- const SAMRAI::hier::Index &fine)
-{
- SAMRAI::hier::Index ip(1,0), jp(0,1);
- SAMRAI::pdat::CellIndex fine_cell(fine);
- SAMRAI::pdat::NodeIndex fine_edge(fine,SAMRAI::pdat::NodeIndex::LowerLeft);
- return (cell(fine_cell) + cell(fine_cell-ip) + cell(fine_cell-jp)
- + cell(fine_cell-ip-jp))/8
- + edge(fine_edge)/4
- + (edge(fine_edge+ip) + edge(fine_edge+jp) + edge(fine_edge-ip)
- + edge(fine_edge-jp))/16;
-}
-
-inline double viscosity_coarsen_3D(const SAMRAI::pdat::CellData<double> &cell,
- const SAMRAI::pdat::EdgeData<double> &edge,
- const SAMRAI::hier::Index &fine)
-{
- SAMRAI::hier::Index ip(1,0,0), jp(0,1,0), kp(0,0,1);
- SAMRAI::pdat::CellIndex fine_cell(fine);
- SAMRAI::pdat::EdgeIndex
- fine_edge_x(fine_cell,SAMRAI::pdat::EdgeIndex::X,
- SAMRAI::pdat::EdgeIndex::LowerLeft),
- fine_edge_y(fine_cell,SAMRAI::pdat::EdgeIndex::Y,
- SAMRAI::pdat::EdgeIndex::LowerLeft),
- fine_edge_z(fine_cell,SAMRAI::pdat::EdgeIndex::Z,
- SAMRAI::pdat::EdgeIndex::LowerLeft);
- return (cell(fine_cell)
- + cell(fine_cell+ip)
- + cell(fine_cell+jp)
- + cell(fine_cell+kp)
- + cell(fine_cell+ip+jp)
- + cell(fine_cell+ip+kp)
- + cell(fine_cell+jp+kp)
- + cell(fine_cell+ip+jp+kp))/32
- + (edge(fine_edge_x+jp+kp) + edge(fine_edge_x+ip+jp+kp)
- + edge(fine_edge_y+ip+kp) + edge(fine_edge_y+ip+jp+kp)
- + edge(fine_edge_z+ip+jp) + edge(fine_edge_z+ip+jp+kp))/8;
-}
-
-#endif
diff -r f089bf615532 -r 567a621add59 wscript
--- a/wscript Thu May 05 13:22:05 2011 -0700
+++ b/wscript Tue Apr 10 11:15:07 2012 -0700
@@ -12,7 +12,7 @@ def build(bld):
features = ['cxx','cprogram'],
source = ['src/main.C',
'src/FACStokes/FACStokes.C',
- 'src/FACStokes/fix_viscosity.C',
+ 'src/FACStokes/fix_moduli.C',
'src/FACStokes/applyGradientDetector.C',
'src/FACStokes/initializeLevelData.C',
'src/FACStokes/packDerivedDataIntoDoubleBuffer.C',
@@ -20,7 +20,6 @@ def build(bld):
'src/FACStokes/setupPlotter.C',
'src/FACStokes/solveStokes.C',
'src/P_Refine.C',
- 'src/P_MDPI_Refine.C',
'src/V_Refine.C',
'src/Resid_Coarsen.C',
'src/V_Coarsen/coarsen_2D.C',
@@ -37,7 +36,6 @@ def build(bld):
'src/StokesFACOps/StokesFACOps.C',
'src/StokesFACOps/computeCompositeResidualOnLevel.C',
'src/StokesFACOps/residual_2D.C',
- 'src/StokesFACOps/residual_3D.C',
'src/StokesFACOps/computeResidualNorm.C',
'src/StokesFACOps/computeVectorWeights.C',
'src/StokesFACOps/deallocateOperatorState.C',
@@ -49,13 +47,9 @@ def build(bld):
'src/StokesFACOps/restrictSolution.C',
'src/StokesFACOps/smoothError.C',
'src/StokesFACOps/smooth_Tackley_2D.C',
- 'src/StokesFACOps/smooth_Tackley_3D.C',
- 'src/StokesFACOps/smooth_Gerya.C',
'src/StokesFACOps/set_boundaries.C',
'src/StokesFACOps/solveCoarsestLevel.C',
- 'src/StokesFACOps/solveCoarsestLevel_HYPRE.C',
'src/StokesFACOps/smooth_V_2D.C',
- 'src/StokesFACOps/smooth_V_3D.C',
'src/StokesFACOps/xeqScheduleGhostFill.C',
'src/StokesFACOps/xeqScheduleGhostFillNoCoarse.C',
'src/StokesFACOps/xeqScheduleProlongation.C',
@@ -77,29 +71,25 @@ def build(bld):
target = 'gamr',
# cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG'],
- cxxflags = ['-g', '-Wall', '-Wextra', '-Wconversion',
+ cxxflags = ['-O3', '-Wall', '-Wextra', '-Wconversion',
'-DTESTING=0'],
lib = ['SAMRAI_appu', 'SAMRAI_algs', 'SAMRAI_solv',
'SAMRAI_geom', 'SAMRAI_mesh', 'SAMRAI_math',
'SAMRAI_pdat', 'SAMRAI_xfer', 'SAMRAI_hier',
- 'SAMRAI_tbox', 'petsc', 'HYPRE',
- 'HYPRE_sstruct_ls', 'HYPRE_sstruct_mv',
- 'HYPRE_struct_ls', 'HYPRE_struct_mv',
- 'HYPRE_parcsr_ls',
- 'HYPRE_DistributedMatrixPilutSolver',
- 'HYPRE_ParaSails', 'HYPRE_Euclid',
- 'HYPRE_MatrixMatrix', 'HYPRE_DistributedMatrix',
- 'HYPRE_IJ_mv', 'HYPRE_parcsr_mv', 'HYPRE_seq_mv',
- 'HYPRE_krylov', 'HYPRE_utilities', 'hdf5',
+ 'SAMRAI_tbox', 'petsc', 'hdf5',
'dl', 'm', 'gfortranbegin', 'gfortran', 'm',
'gfortranbegin', 'gfortran', 'm', 'gfortranbegin',
- 'gfortran', 'm', 'sundials_cvode', 'sundials_kinsol'],
- libpath = ['../../SAMRAI-hg-build/lib',
- '/usr/lib/petsc/linux-gnu-c-opt/lib'],
- includes = ['src','../SAMRAI-hg-build/include',
- '../SAMRAI-hg/source',
- '/usr/lib/petsc/include',
+ 'gfortran', 'm'],
+ libpath = ['/Users/sbarbot/Documents/src/samrai/lib',
+ '/Users/sbarbot/Documents/src/petsc/petsc-3.2-p7/real8-with-debug/lib',
+ '/sw/lib',
+ '/sw/lib/gcc4.6/lib/gcc/x86_64-apple-darwin10.8.0/4.6.1',
+ '/sw/lib/gcc4.6/lib',
+ '/shared/lib'],
+ includes = ['src','../../samrai/include',
+ '../../samrai/source',
+ '../../petsc/petsc-3.2-p7/real8-with-debug/include',
'/usr/lib/petsc/bmake/linux-gnu-c-opt',
- '/usr/lib/petsc/src/vec',
+ '../../petsc/petsc-3.2-p7/src/vec',
'/usr/include']
)
More information about the CIG-COMMITS
mailing list