[cig-commits] commit: Change adaptation to adaption internally. Add the ability to specify v_rhs in the input file. Add the ability to specify minimum level of total refinement.
Mercurial
hg at geodynamics.org
Mon May 2 12:25:16 PDT 2011
changeset: 227:5864f64ab798
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Mon May 02 12:23:53 2011 -0700
files: src/FACStokes.h src/FACStokes/FACStokes.C src/FACStokes/applyGradientDetector.C src/FACStokes/initializeLevelData.C
description:
Change adaptation to adaption internally. Add the ability to specify v_rhs in the input file. Add the ability to specify minimum level of total refinement.
diff -r 08ca3b14e5fc -r 5864f64ab798 src/FACStokes.h
--- a/src/FACStokes.h Mon May 02 10:03:56 2011 -0700
+++ b/src/FACStokes.h Mon May 02 12:23:53 2011 -0700
@@ -182,14 +182,17 @@ namespace SAMRAI {
* These are initialized in the constructor and never change.
*/
- double d_adaptation_threshold;
+ double d_adaption_threshold;
+ int min_full_refinement_level;
public:
int p_id, cell_viscosity_id, edge_viscosity_id, dp_id, p_exact_id,
p_rhs_id, v_id, v_rhs_id;
- tbox::Array<double> viscosity;
+ tbox::Array<double> viscosity, viscosity_xyz_max, viscosity_xyz_min;
tbox::Array<int> viscosity_ijk;
- tbox::Array<double> viscosity_xyz_max, viscosity_xyz_min;
+
+ tbox::Array<double> v_rhs, v_rhs_xyz_max, v_rhs_xyz_min;
+ tbox::Array<int> v_rhs_ijk;
//@}
};
diff -r 08ca3b14e5fc -r 5864f64ab798 src/FACStokes/FACStokes.C
--- a/src/FACStokes/FACStokes.C Mon May 02 10:03:56 2011 -0700
+++ b/src/FACStokes/FACStokes.C Mon May 02 12:23:53 2011 -0700
@@ -63,16 +63,17 @@ namespace SAMRAI {
*/
tbox::Pointer<pdat::CellVariable<double> >
- p(new pdat::CellVariable<double>(dim, object_name + ":p", 1));
- p_id = vdb->registerVariableAndContext(p, d_context, hier::IntVector(dim, 1)
+ p_ptr(new pdat::CellVariable<double>(dim, object_name + ":p", 1));
+ p_id = vdb->registerVariableAndContext(p_ptr, d_context,
+ hier::IntVector(dim, 1)
/* ghost cell width is 1 for
stencil widths */);
tbox::Pointer<pdat::CellVariable<double> >
- cell_viscosity(new pdat::CellVariable<double>(dim,
- object_name
- + ":cell_viscosity"));
- cell_viscosity_id = vdb->registerVariableAndContext(cell_viscosity,
+ cell_viscosity_ptr(new pdat::CellVariable<double>(dim,
+ object_name
+ + ":cell_viscosity"));
+ cell_viscosity_id = vdb->registerVariableAndContext(cell_viscosity_ptr,
d_context,
hier::IntVector(dim, 1)
/* ghost cell width is
@@ -81,11 +82,11 @@ namespace SAMRAI {
if(dim.getValue()==2)
{
tbox::Pointer<pdat::NodeVariable<double> >
- edge_viscosity(new pdat::NodeVariable<double>(dim,
- object_name
- + ":edge_viscosity"));
+ edge_viscosity_ptr(new pdat::NodeVariable<double>(dim,
+ object_name
+ + ":edge_viscosity"));
edge_viscosity_id =
- vdb->registerVariableAndContext(edge_viscosity,d_context,
+ vdb->registerVariableAndContext(edge_viscosity_ptr,d_context,
hier::IntVector(dim,1)
/* ghost cell width is 1 in
case needed */);
@@ -93,53 +94,56 @@ namespace SAMRAI {
else if(dim.getValue()==3)
{
tbox::Pointer<pdat::EdgeVariable<double> >
- edge_viscosity(new pdat::EdgeVariable<double>(dim,
- object_name
- + ":edge_viscosity"));
+ edge_viscosity_ptr(new pdat::EdgeVariable<double>(dim,
+ object_name
+ + ":edge_viscosity"));
edge_viscosity_id =
- vdb->registerVariableAndContext(edge_viscosity,d_context,
+ vdb->registerVariableAndContext(edge_viscosity_ptr,d_context,
hier::IntVector(dim,1)
/* ghost cell width is 1 in
case needed */);
}
tbox::Pointer<pdat::CellVariable<double> >
- dp(new pdat::CellVariable<double>(dim, object_name + ":dp"));
- dp_id = vdb->registerVariableAndContext(dp,d_context,
+ dp_ptr(new pdat::CellVariable<double>(dim, object_name + ":dp"));
+ dp_id = vdb->registerVariableAndContext(dp_ptr,d_context,
hier::IntVector(dim, 1)
/* ghost cell width is
1 in case needed */);
tbox::Pointer<pdat::CellVariable<double> >
- p_exact(new pdat::CellVariable<double>(dim, object_name + ":p exact"));
- p_exact_id = vdb->registerVariableAndContext(p_exact,d_context,
+ p_exact_ptr(new pdat::CellVariable<double>(dim, object_name + ":p exact"));
+ p_exact_id = vdb->registerVariableAndContext(p_exact_ptr,d_context,
hier::IntVector(dim, 1)
/* ghost cell width is
1 in case needed */);
tbox::Pointer<pdat::CellVariable<double> >
- p_rhs(new pdat::CellVariable<double>(dim,object_name
- + ":p right hand side"));
- p_rhs_id = vdb->registerVariableAndContext(p_rhs,d_context,
+ p_rhs_ptr(new pdat::CellVariable<double>(dim,object_name
+ + ":p right hand side"));
+ p_rhs_id = vdb->registerVariableAndContext(p_rhs_ptr,d_context,
hier::IntVector(dim, 1));
tbox::Pointer<pdat::SideVariable<double> >
- v(new pdat::SideVariable<double>(dim, object_name + ":v", 1));
- v_id = vdb->registerVariableAndContext(v, d_context, hier::IntVector(dim, 1)
+ v_ptr(new pdat::SideVariable<double>(dim, object_name + ":v", 1));
+ v_id = vdb->registerVariableAndContext(v_ptr, d_context,
+ hier::IntVector(dim, 1)
/* ghost cell width is 1 for
stencil widths */);
tbox::Pointer<pdat::SideVariable<double> >
- v_rhs(new pdat::SideVariable<double>(dim,object_name
- + ":v right hand side"));
- v_rhs_id = vdb->registerVariableAndContext(v_rhs,d_context,
+ v_rhs_ptr(new pdat::SideVariable<double>(dim,object_name
+ + ":v right hand side"));
+ v_rhs_id = vdb->registerVariableAndContext(v_rhs_ptr,d_context,
hier::IntVector(dim, 1)
/* ghost cell width is
1 for coarsening
operator */);
- d_adaptation_threshold=database->getDoubleWithDefault("adaption_threshold",
- 1.0e-15);
+ d_adaption_threshold=database->getDoubleWithDefault("adaption_threshold",
+ 1.0e-15);
+ min_full_refinement_level
+ =database->getIntegerWithDefault("min_full_refinement_level",0);
if(database->keyExists("viscosity_data"))
{
@@ -147,6 +151,14 @@ namespace SAMRAI {
viscosity_xyz_min=database->getDoubleArray("viscosity_coord_min");
viscosity_xyz_max=database->getDoubleArray("viscosity_coord_max");
viscosity=database->getDoubleArray("viscosity_data");
+ }
+
+ if(database->keyExists("v_rhs_data"))
+ {
+ v_rhs_ijk=database->getIntegerArray("v_rhs_ijk");
+ v_rhs_xyz_min=database->getDoubleArray("v_rhs_coord_min");
+ v_rhs_xyz_max=database->getDoubleArray("v_rhs_coord_max");
+ v_rhs=database->getDoubleArray("v_rhs_data");
}
/*
diff -r 08ca3b14e5fc -r 5864f64ab798 src/FACStokes/applyGradientDetector.C
--- a/src/FACStokes/applyGradientDetector.C Mon May 02 10:03:56 2011 -0700
+++ b/src/FACStokes/applyGradientDetector.C Mon May 02 12:23:53 2011 -0700
@@ -65,14 +65,15 @@ void SAMRAI::FACStokes::applyGradientDet
// << (estimate_data(cell_index) > d_adaptation_threshold)
// << " "
// << "\n";
- if (estimate_data(cell_index) > d_adaptation_threshold)
+ if (estimate_data(cell_index) > d_adaption_threshold
+ || ln<min_full_refinement_level)
{
tag_cell_data(cell_index) = 1;
++ntag;
}
}
}
- tbox::plog << "Adaption threshold is " << d_adaptation_threshold << "\n";
+ tbox::plog << "Adaption threshold is " << d_adaption_threshold << "\n";
tbox::plog << "Number of cells tagged on level " << ln << " is "
<< ntag << "/" << ntotal << "\n";
tbox::plog << "Max estimate is " << maxestimate << "\n";
diff -r 08ca3b14e5fc -r 5864f64ab798 src/FACStokes/initializeLevelData.C
--- a/src/FACStokes/initializeLevelData.C Mon May 02 10:03:56 2011 -0700
+++ b/src/FACStokes/initializeLevelData.C Mon May 02 12:23:53 2011 -0700
@@ -33,6 +33,7 @@ void SAMRAI::FACStokes::initializeLevelD
tbox::Pointer<hier::PatchLevel> level =
hierarchy->getPatchLevel(level_number);
+ const int dim=d_dim.getValue();
if (allocate_data) {
level->allocatePatchData(p_id);
@@ -44,13 +45,6 @@ void SAMRAI::FACStokes::initializeLevelD
level->allocatePatchData(v_id);
level->allocatePatchData(v_rhs_id);
}
-
- // const double inclusion_radius=0.5;
- const double inclusion_viscosity=1e2;
- const double background_viscosity=1;
-
- // const double background_density(1), block_density(1);
- const double background_density(1), block_density(1.03);
/*
* Initialize data in all patches in the level.
@@ -65,7 +59,9 @@ void SAMRAI::FACStokes::initializeLevelD
}
tbox::Pointer<geom::CartesianPatchGeometry>
geom = patch->getPatchGeometry();
+ const double *dx=geom->getDx();
+ /* Initialize cell viscosity */
tbox::Pointer<pdat::CellData<double> > cell_viscosity_ptr =
patch->getPatchData(cell_viscosity_id);
pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
@@ -79,110 +75,80 @@ void SAMRAI::FACStokes::initializeLevelD
double y=geom->getXLower()[1]
+ geom->getDx()[1]*(c[1]-cell_visc_box.lower()[1] + 0.5);
- if(!viscosity.empty())
+ int i(static_cast<int>(x*(viscosity_ijk[0]-1)
+ /(viscosity_xyz_max[0]-viscosity_xyz_min[0]))),
+ j(static_cast<int>(y*(viscosity_ijk[1]-1)
+ /(viscosity_xyz_max[1]-viscosity_xyz_min[1])));
+ i=std::max(0,std::min(viscosity_ijk[0]-1,i));
+ j=std::max(0,std::min(viscosity_ijk[1]-1,j));
+
+ if(dim==2)
{
- int i(static_cast<int>(x*(viscosity_ijk[0]-1)
- /(viscosity_xyz_max[0]-viscosity_xyz_min[0]))),
- j(static_cast<int>(y*(viscosity_ijk[1]-1)
- /(viscosity_xyz_max[1]-viscosity_xyz_min[1])));
- i=std::max(0,std::min(viscosity_ijk[0]-1,i));
- j=std::max(0,std::min(viscosity_ijk[1]-1,j));
-
- if(d_dim.getValue()==2)
- cell_viscosity(c)=viscosity[i+viscosity_ijk[0]*j];
- else
- {
- double z=geom->getXLower()[2]
- + geom->getDx()[2]*(c[2]-cell_visc_box.lower()[2] + 0.5);
- int k(static_cast<int>(z*(viscosity_ijk[2]-1)
- /(viscosity_xyz_max[2]
- - viscosity_xyz_min[2])));
- k=std::max(0,std::min(viscosity_ijk[2]-1,k));
- cell_viscosity(c)=
- viscosity[i+viscosity_ijk[0]*(j+viscosity_ijk[1]*k)];
- }
+ cell_viscosity(c)=viscosity[i+viscosity_ijk[0]*j];
}
else
{
- // if(x*x + y*y < inclusion_radius*inclusion_radius)
- // cell_viscosity(c)=inclusion_viscosity;
- // else
- // cell_viscosity(c)=background_viscosity;
-
- if(d_dim.getValue()==2)
- {
- if(x<1.0/3 || x>2.0/3 || y<1.0/3 || y>2.0/3)
- cell_viscosity(c)=background_viscosity;
- else
- cell_viscosity(c)=inclusion_viscosity;
- }
- else
- {
- double z=geom->getXLower()[2]
- + geom->getDx()[2]*(c[2]-cell_visc_box.lower()[2] + 0.5);
- if(x<1.0/3 || x>2.0/3 || y<1.0/3 || y>2.0/3 || z<1.0/3 || z>2.0/3)
- cell_viscosity(c)=background_viscosity;
- else
- cell_viscosity(c)=inclusion_viscosity;
- }
+ double z=geom->getXLower()[2]
+ + geom->getDx()[2]*(c[2]-cell_visc_box.lower()[2] + 0.5);
+ int k(static_cast<int>(z*(viscosity_ijk[2]-1)
+ /(viscosity_xyz_max[2]
+ - viscosity_xyz_min[2])));
+ k=std::max(0,std::min(viscosity_ijk[2]-1,k));
+ cell_viscosity(c)=
+ viscosity[i+viscosity_ijk[0]*(j+viscosity_ijk[1]*k)];
}
}
+ /* I do not think this is actually necessary. */
tbox::Pointer<pdat::CellData<double> > dp_data =
patch->getPatchData(dp_id);
-
- /* I do not think this is actually necessary. */
dp_data->fill(0.0);
tbox::Pointer<pdat::CellData<double> > p_rhs_data =
patch->getPatchData(p_rhs_id);
-
p_rhs_data->fill(0.0);
+ /* v_rhs */
tbox::Pointer<pdat::SideData<double> > v_rhs_data =
patch->getPatchData(v_rhs_id);
- hier::Box pbox = v_rhs_data->getBox();
+ if(v_rhs.empty())
+ {
+ v_rhs_data->fill(0,0);
+ }
+ else
+ {
+ hier::Box pbox = v_rhs_data->getBox();
+ int ix_offset(0);
+ for(int ix=0;ix<dim;++ix)
+ {
+ double offset[]={0,0,0};
+ offset[ix]=0.5;
- for(pdat::SideIterator si(pbox,0); si; si++)
- {
- pdat::SideIndex s=si();
- (*v_rhs_data)(s)=0;
+ for(pdat::SideIterator si(pbox,ix); si; si++)
+ {
+ pdat::SideIndex s=si();
+ double xyz[]={0,0,0};
+ for(int d=0;d<dim;++d)
+ xyz[d]=geom->getXLower()[d]
+ + dx[d]*(s[d]-pbox.lower()[d]+offset[d]);
+
+ int ijk(0), factor(1);
+ for(int d=0;d<dim;++d)
+ {
+ int i=static_cast<int>(xyz[d]*(v_rhs_ijk[d]-1)
+ /(v_rhs_xyz_max[d]-v_rhs_xyz_min[d]));
+ i=std::max(0,std::min(v_rhs_ijk[d]-1,i));
+ ijk+=i*factor;
+ factor*=v_rhs_ijk[d];
+ }
+ (*v_rhs_data)(s)=v_rhs[ijk+ix_offset];
+ }
+ int i=1;
+ for(int d=0;d<dim;++d)
+ i*=v_rhs_ijk[d];
+ ix_offset+=i;
+ }
}
- for(pdat::SideIterator si(pbox,1); si; si++)
- {
- pdat::SideIndex s=si();
-
- double x=geom->getXLower()[0]
- + geom->getDx()[0]*(s[0]-pbox.lower()[0]+0.5);
- double y=geom->getXLower()[1]
- + geom->getDx()[1]*(s[1]-pbox.lower()[1]);
-
- if(d_dim.getValue()==2)
- {
- if(x<1.0/3 || x>2.0/3 || y<1.0/3 || y>2.0/3)
- (*v_rhs_data)(s)=background_density;
- else
- (*v_rhs_data)(s)=block_density;
- }
- else
- {
- double z=geom->getXLower()[2]
- + geom->getDx()[2]*(s[2]-pbox.lower()[2]);
- if(x<1.0/3 || x>2.0/3 || y<1.0/3 || y>2.0/3 || z<1.0/3 || z>2.0/3)
- (*v_rhs_data)(s)=background_density;
- else
- (*v_rhs_data)(s)=block_density;
- }
-
- // (*v_rhs_data)(s)=10;
- // (*v_rhs_data)(s)=0;
- }
- if(d_dim.getValue()==3)
- for(pdat::SideIterator si(pbox,2); si; si++)
- {
- pdat::SideIndex s=si();
- (*v_rhs_data)(s)=0;
- }
} // End patch loop.
}
More information about the CIG-COMMITS
mailing list