[cig-commits] commit: 3D on a single level seems to be working in parallel
Mercurial
hg at geodynamics.org
Sun Mar 20 00:11:30 PDT 2011
changeset: 145:6e0fdd8950bd
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sun Mar 20 00:10:31 2011 -0700
files: Edge_Viscosity_Coarsen.h FACStokes/FACStokes.C FACStokes/initializeLevelData.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/residual_3D.C input/shear_corner_3D.input
description:
3D on a single level seems to be working in parallel
diff -r cb8e50485149 -r 6e0fdd8950bd Edge_Viscosity_Coarsen.h
--- a/Edge_Viscosity_Coarsen.h Sat Mar 19 20:15:50 2011 -0700
+++ b/Edge_Viscosity_Coarsen.h Sun Mar 20 00:10:31 2011 -0700
@@ -20,6 +20,7 @@
#include "SAMRAI/hier/Patch.h"
#include "SAMRAI/tbox/Pointer.h"
#include "SAMRAI/pdat/NodeVariable.h"
+#include "SAMRAI/pdat/EdgeVariable.h"
#include <string>
@@ -66,8 +67,9 @@ public:
{
TBOX_DIM_ASSERT_CHECK_ARGS2(*this, *var);
- const tbox::Pointer<pdat::NodeVariable<double> > cast_var(var);
- if (!cast_var.isNull() && (op_name == d_name_id)) {
+ const tbox::Pointer<pdat::NodeVariable<double> > node_var(var);
+ const tbox::Pointer<pdat::EdgeVariable<double> > edge_var(var);
+ if ((!node_var.isNull() || !edge_var.isNull()) && (op_name == d_name_id)) {
return true;
} else {
return false;
diff -r cb8e50485149 -r 6e0fdd8950bd FACStokes/FACStokes.C
--- a/FACStokes/FACStokes.C Sat Mar 19 20:15:50 2011 -0700
+++ b/FACStokes/FACStokes.C Sun Mar 20 00:10:31 2011 -0700
@@ -79,16 +79,30 @@ namespace SAMRAI {
/* ghost cell width is
1 in case needed */);
- tbox::Pointer<pdat::NodeVariable<double> >
- edge_viscosity(new pdat::NodeVariable<double>(dim,
- object_name
- + ":edge_viscosity"));
- edge_viscosity_id = vdb->registerVariableAndContext(edge_viscosity,
- d_context,
- hier::IntVector(dim, 1)
- /* ghost cell width is
- 1 in case needed */);
-
+ if(dim.getValue()==2)
+ {
+ tbox::Pointer<pdat::NodeVariable<double> >
+ edge_viscosity(new pdat::NodeVariable<double>(dim,
+ object_name
+ + ":edge_viscosity"));
+ edge_viscosity_id =
+ vdb->registerVariableAndContext(edge_viscosity,d_context,
+ hier::IntVector(dim,1)
+ /* ghost cell width is 1 in
+ case needed */);
+ }
+ else if(dim.getValue()==3)
+ {
+ tbox::Pointer<pdat::EdgeVariable<double> >
+ edge_viscosity(new pdat::EdgeVariable<double>(dim,
+ object_name
+ + ":edge_viscosity"));
+ edge_viscosity_id =
+ vdb->registerVariableAndContext(edge_viscosity,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"));
diff -r cb8e50485149 -r 6e0fdd8950bd FACStokes/initializeLevelData.C
--- a/FACStokes/initializeLevelData.C Sat Mar 19 20:15:50 2011 -0700
+++ b/FACStokes/initializeLevelData.C Sun Mar 20 00:10:31 2011 -0700
@@ -84,83 +84,68 @@ void SAMRAI::FACStokes::initializeLevelD
patch->getPatchData(cell_viscosity_id);
pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
- hier::Box visc_box = cell_viscosity.getBox();
+ hier::Box cell_visc_box = cell_viscosity.getBox();
for(pdat::CellIterator ci(cell_viscosity.getGhostBox()); ci; ci++)
{
pdat::CellIndex c=ci();
double x=geom->getXLower()[0]
- + geom->getDx()[0]*(c[0]-visc_box.lower()[0] + 0.5);
+ + geom->getDx()[0]*(c[0]-cell_visc_box.lower()[0] + 0.5);
double y=geom->getXLower()[1]
- + geom->getDx()[1]*(c[1]-visc_box.lower()[1] + 0.5);
+ + geom->getDx()[1]*(c[1]-cell_visc_box.lower()[1] + 0.5);
if(x*x + y*y < inclusion_radius*inclusion_radius)
cell_viscosity(c)=inclusion_viscosity;
else
cell_viscosity(c)=background_viscosity;
-
- // tbox::plog << "cell "
- // << c[0] << " "
- // << c[1] << " "
- // << x << " "
- // << y << " "
- // << geom->getXLower()[0] << " "
- // << geom->getXLower()[1] << " "
- // << geom->getDx()[0] << " "
- // << geom->getDx()[1] << " "
- // << std::boolalpha
- // << (x*x + y*y < inclusion_radius*inclusion_radius) << " "
- // << cell_viscosity(c) << " "
- // << "\n";
-
}
- tbox::Pointer<pdat::NodeData<double> > edge_viscosity_ptr =
- patch->getPatchData(edge_viscosity_id);
- pdat::NodeData<double> &edge_viscosity(*edge_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);
+ hier::Box edge_visc_box = edge_viscosity.getBox();
- for(pdat::NodeIterator ei(edge_viscosity.getGhostBox()); ei; ei++)
+ for(pdat::NodeIterator ei(edge_viscosity.getGhostBox()); ei; ei++)
+ {
+ pdat::NodeIndex e=ei();
+ double x=geom->getXLower()[0]
+ + geom->getDx()[0]*(e[0]-edge_visc_box.lower()[0]);
+ double y=geom->getXLower()[1]
+ + geom->getDx()[1]*(e[1]-edge_visc_box.lower()[1]);
+ if(x*x + y*y < inclusion_radius*inclusion_radius)
+ edge_viscosity(e)=inclusion_viscosity;
+ else
+ edge_viscosity(e)=background_viscosity;
+ }
+ }
+ else if(d_dim.getValue()==3)
{
- pdat::NodeIndex e=ei();
- double x=geom->getXLower()[0]
- + geom->getDx()[0]*(e[0]-visc_box.lower()[0]);
- double y=geom->getXLower()[1]
- + geom->getDx()[1]*(e[1]-visc_box.lower()[1]);
- if(x*x + y*y < inclusion_radius*inclusion_radius)
- // if(x<inclusion_radius && y<inclusion_radius)
- edge_viscosity(e)=inclusion_viscosity;
- else
- edge_viscosity(e)=background_viscosity;
+ tbox::Pointer<pdat::EdgeData<double> > edge_viscosity_ptr =
+ patch->getPatchData(edge_viscosity_id);
+ pdat::EdgeData<double> &edge_viscosity(*edge_viscosity_ptr);
+ hier::Box edge_visc_box = edge_viscosity.getBox();
- // tbox::plog << "edge "
- // << e[0] << " "
- // << e[1] << " "
- // << x << " "
- // << y << " "
- // << geom->getXLower()[0] << " "
- // << geom->getXLower()[1] << " "
- // << geom->getDx()[0] << " "
- // << geom->getDx()[1] << " "
- // << std::boolalpha
- // << (x*x + y*y < inclusion_radius*inclusion_radius) << " "
- // << edge_viscosity(e) << " "
- // << "\n";
+ for(int ix=0;ix<3;++ix)
+ for(pdat::EdgeIterator ei(edge_viscosity.getGhostBox(),ix); ei; ei++)
+ {
+ pdat::EdgeIndex e=ei();
+ double dx(0);
+ if(ix==0)
+ dx=0.5;
+ double dy(0);
+ if(ix==1)
+ dy=0.5;
+ double x=geom->getXLower()[0]
+ + geom->getDx()[0]*(e[0]-edge_visc_box.lower()[0] + dx);
+ double y=geom->getXLower()[1]
+ + geom->getDx()[1]*(e[1]-edge_visc_box.lower()[1] + dy);
+ if(x*x + y*y < inclusion_radius*inclusion_radius)
+ edge_viscosity(e)=inclusion_viscosity;
+ else
+ edge_viscosity(e)=background_viscosity;
+ }
}
-
-
- // tbox::Pointer<pdat::CellData<double> > cell_viscosity_data =
- // patch->getPatchData(cell_viscosity_id);
-
- // /* At some point this needs to do the proper interpolation for
- // lower levels */
- // cell_viscosity_data->fill(1.0);
-
- // tbox::Pointer<pdat::NodeData<double> > edge_viscosity_data =
- // patch->getPatchData(edge_viscosity_id);
-
- // /* At some point this needs to do the proper interpolation for
- // lower levels */
- // edge_viscosity_data->fill(1.0);
-
tbox::Pointer<pdat::CellData<double> > dp_data =
patch->getPatchData(dp_id);
@@ -186,8 +171,14 @@ void SAMRAI::FACStokes::initializeLevelD
for(pdat::SideIterator si(pbox,1); si; si++)
{
pdat::SideIndex s=si();
- // (*v_rhs_data)(s)=10;
- (*v_rhs_data)(s)=0;
+ (*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.
}
diff -r cb8e50485149 -r 6e0fdd8950bd StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Sat Mar 19 20:15:50 2011 -0700
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Sun Mar 20 00:10:31 2011 -0700
@@ -95,7 +95,8 @@ void SAMRAI::solv::StokesFACOps::compute
*p_resid_ptr,*v_resid_ptr,*patch,pbox,*geom);
break;
case 3:
- // residual_3D();
+ residual_3D(*p_ptr,*v_ptr,*cell_viscosity_ptr,*p_rhs_ptr,*v_rhs_ptr,
+ *p_resid_ptr,*v_resid_ptr,*patch,pbox,*geom);
break;
default:
abort();
diff -r cb8e50485149 -r 6e0fdd8950bd StokesFACOps/residual_3D.C
--- a/StokesFACOps/residual_3D.C Sat Mar 19 20:15:50 2011 -0700
+++ b/StokesFACOps/residual_3D.C Sun Mar 20 00:10:31 2011 -0700
@@ -13,15 +13,13 @@ void SAMRAI::solv::StokesFACOps::residua
const hier::Box &pbox,
const geom::CartesianPatchGeometry &geom)
{
- const hier::Index ip(1,0,0), jp(0,1,0), kp(0,0,1);
-
tbox::Pointer<pdat::EdgeData<double> >
edge_viscosity_ptr = patch.getPatchData(edge_viscosity_id);
pdat::EdgeData<double> &edge_viscosity(*edge_viscosity_ptr);
- double dx = geom.getDx()[0];
- double dy = geom.getDx()[1];
- double dz = geom.getDx()[2];
+ const double *Dx = geom.getDx();
+ const hier::Index ip(1,0,0), jp(0,1,0), kp(0,0,1);
+ const hier::Index pp[]={ip,jp,kp};
for(pdat::CellIterator ci(pbox); ci; ci++)
{
@@ -36,75 +34,48 @@ void SAMRAI::solv::StokesFACOps::residua
++front[2];
--back[2];
- const pdat::SideIndex
- x(center,0,pdat::SideIndex::Lower),
- y(center,1,pdat::SideIndex::Lower),
- z(center,2,pdat::SideIndex::Lower);
- const pdat::EdgeIndex
- edge_x(center,0,pdat::EdgeIndex::LowerLeft),
- edge_y(center,1,pdat::EdgeIndex::LowerLeft),
- edge_z(center,2,pdat::EdgeIndex::LowerLeft);
-
/* p */
if(center[0]!=pbox.upper(0) && center[1]!=pbox.upper(1)
&& center[2]!=pbox.upper(2))
{
- double dvx_dx=(v(x+ip) - v(x))/dx;
- double dvy_dy=(v(y+jp) - v(y))/dy;
- double dvz_dz=(v(z+kp) - v(z))/dy;
+ 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;
}
- /* vx */
- if(center[1]!=pbox.upper(1) && center[2]!=pbox.upper(2))
+ for(int ix=0;ix<3;++ix)
{
- /* If x==0 */
- if((center[0]==pbox.lower(0) && v(x-ip)==boundary_value)
- || (center[0]==pbox.upper(0) && v(x+ip)==boundary_value))
+ const int iy((ix+1)%3), iz((ix+2)%3);
+ const pdat::SideIndex
+ x(center,ix,pdat::SideIndex::Lower),
+ y(center,iy,pdat::SideIndex::Lower),
+ z(center,iz,pdat::SideIndex::Lower);
+ const pdat::EdgeIndex
+ edge_y(center,iy,pdat::EdgeIndex::LowerLeft),
+ edge_z(center,iz,pdat::EdgeIndex::LowerLeft);
+
+ if(center[iy]!=pbox.upper(iy) && center[iz]!=pbox.upper(iz))
{
- v_resid(x)=0;
+ if((center[ix]==pbox.lower(ix) && v(x-pp[ix])==boundary_value)
+ || (center[ix]==pbox.upper(ix) && v(x+pp[ix])==boundary_value))
+ {
+ v_resid(x)=0;
+ }
+ else
+ {
+ v_resid(x)=v_rhs(x)
+ - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
+ center,edge_y,edge_z,x,y,z,
+ pp[ix],pp[iy],pp[iz],Dx[ix],Dx[iy],Dx[iz]);
+ }
}
- else
- {
- v_resid(x)=v_rhs(x)
- - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
- center,edge_y,edge_z,x,y,z,ip,jp,kp,dx,dy,dz);
- }
- }
-
- /* vy */
- if(center[0]!=pbox.upper(0) && center[2]!=pbox.upper(2))
- {
- /* If y==0 */
- if((center[1]==pbox.lower(1) && v(y-jp)==boundary_value)
- || (center[1]==pbox.upper(1) && v(y+jp)==boundary_value))
- {
- v_resid(y)=0;
- }
- else
- {
- v_resid(y)=v_rhs(y)
- - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
- center,edge_z,edge_x,y,z,x,jp,kp,ip,dy,dz,dx);
- }
- }
-
- /* vz */
- if(center[1]!=pbox.upper(0) && center[2]!=pbox.upper(2))
- {
- /* If z==0 */
- if((center[2]==pbox.lower(1) && v(z-kp)==boundary_value)
- || (center[2]==pbox.upper(1) && v(z+kp)==boundary_value))
- {
- v_resid(z)=0;
- }
- else
- {
- v_resid(z)=v_rhs(z)
- - v_operator_3D(v,p,cell_viscosity,edge_viscosity,
- center,edge_x,edge_y,z,x,y,kp,ip,jp,dz,dx,dy);
- }
- }
+ }
}
}
diff -r cb8e50485149 -r 6e0fdd8950bd input/shear_corner_3D.input
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/input/shear_corner_3D.input Sun Mar 20 00:10:31 2011 -0700
@@ -0,0 +1,137 @@
+/*************************************************************************
+ *
+ * This file is part of the SAMRAI distribution. For full copyright
+ * information, see COPYRIGHT and COPYING.LESSER.
+ *
+ * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
+ * Description: Input file for example FAC Stokes solver
+ *
+ ************************************************************************/
+// This is the input file for the 2D FAC example
+// demonstrating changes in boundary conditions.
+//
+// Note that using constant refine prolongation
+// lead to slower convergence.
+
+Main {
+ dim = 3
+ // Base name for output files.
+ base_name = "shear_corner"
+ // Whether to log all nodes in a parallel run.
+ log_all_nodes = TRUE
+ // Visualization writers to write files for.
+ vis_writer = "Vizamrai", "VisIt"
+}
+
+FACStokes {
+ // The FACStokes class is the "user class" in this example.
+ // It owns the solver and contains the code to set up the solver.
+ // The inputs for FACStokes is simply the inputs for the individual
+ // parts owned by the FACStokes class.
+ fac_solver {
+ // This is the input for the cell-centered Stokes FAC solver
+ // class in the SAMRAI library.
+ enable_logging = TRUE // Bool flag to switch logging on/off
+ max_cycles = 100 // Max number of FAC cycles to use
+ residual_tol = 1e-8 // Residual tolerance to solve for
+ num_pre_sweeps = 5 // Number of presmoothing sweeps to use
+ num_post_sweeps = 5 // Number of postsmoothing sweeps to use
+ smoothing_choice = "Tackley"
+ coarse_solver_choice = "Tackley"
+ // smoothing_choice = "Gerya"
+ // coarse_solver_choice = "Gerya"
+ coarse_solver_max_iterations = 25
+ coarse_solver_tolerance = 1e-8
+ // p_prolongation_method = "P_MDPI_REFINE"
+ p_prolongation_method = "P_REFINE"
+ v_prolongation_method = "V_REFINE"
+ }
+ bc_coefs {
+ // These are the boundary condition specifications. The number
+ // after "boundary_" is the location index of the boundary.
+ // The inputs are arrays of strings where the first string
+ // indicates the type of values you want to set. "slope" means
+ // boundary slope, "value" means boundary value, and "coefficients"
+ // mean the raw Robin boundary condition coefficients.
+ // The remaining strings are converted into numbers as
+ // appropriate for what boundary condition you specified with
+ // the first string. Other boundary conditions are possible.
+ // see the solv_RobinBcCoefStrategy class.
+ // Examples:
+ // boundary_0 = "slope", "0"
+ // boundary_1 = "coefficients", "0", "0"
+ // boundary_2 = "value", "0"
+ // boundary_3 = "value", "0"
+ boundary_0 = "value", "0"
+ boundary_1 = "value", "0"
+ boundary_2 = "value", "0"
+ boundary_3 = "value", "0"
+ }
+}
+
+CartesianGridGeometry {
+ // Specify lower/upper corners of the computational domain and a
+ // set of non-overlapping boxes defining domain interior. If union
+ // of boxes is not a parallelpiped, lower/upper corner data corresponds
+ // to min/max corner indices over all boxes given.
+ // x_lo -- (double array) lower corner of computational domain [REQD]
+ // x_up -- (double array) upper corner of computational domain [REQD]
+ // domain_boxes -- (box array) set of boxes that define interior of
+ // physical domain. [REQD]
+ // periodic_dimension -- (int array) coordinate directions in which
+ // domain is periodic. Zero indicates not
+ // periodic, non-zero value indicates periodicity.
+ // [0]
+ domain_boxes = [(0,0,0), (7,7,7)]
+ x_lo = 0, 0, 0
+ x_up = 1, 1, 1
+}
+
+StandardTagAndInitialize {
+ tagging_method = "REFINE_BOXES"
+ RefineBoxes {
+ level_0 = [(0,0,0),(3,3,3)]
+ }
+}
+
+PatchHierarchy {
+ // Information used to create patches in AMR hierarchy.
+ // max_levels -- (int) max number of mesh levels in hierarchy [REQD]
+ //
+ // For most of the following parameters, the number of precribed data
+ // values need not match the number of levels in the hierarchy
+ // (determined by max_levels). If more values are given than number
+ // of levels, extraneous values will be ignored. If less are give, then
+ // values that correspond to individual levels will apply to those
+ // levels. Missing values will be taken from those for the finest
+ // level specified.
+ //
+ // ratio_to_coarser {
+ // level_1 -- (int array) ratio between index spaces on
+ // level 1 to level 0 [REQD]
+ // level_2 -- (int array) ratio between index spaces on
+ // level 2 to level 1 [REQD]
+ // etc....
+ // }
+ // largest_patch_size {
+ // level_0 -- (int array) largest patch allowed on level 0.
+ // [REQD]
+ // level_1 -- (int array) " " " " level 1
+ // [level 0 entry]
+ // etc....
+ // }
+ max_levels = 1
+ ratio_to_coarser {
+ level_1 = 2, 2, 2
+ level_2 = 2, 2, 2
+ level_3 = 2, 2, 2
+ level_4 = 2, 2, 2
+ }
+ largest_patch_size {
+ level_0 = 8, 8, 8
+ // all finer levels will use same values as level_0...
+ }
+}
+
+GriddingAlgorithm {
+}
More information about the CIG-COMMITS
mailing list