[cig-commits] commit: Experiment with new coarsening and refinement operators. In progress and does not really work.
Mercurial
hg at geodynamics.org
Thu Jun 7 13:36:02 PDT 2012
changeset: 276:b6526d45fde8
user: Walter Landry <wlandry at caltech.edu>
date: Tue Jun 05 14:15:49 2012 -0700
files: src/V_Coarsen/coarsen_2D.C src/V_Refine.C src/V_Refine.h src/V_Refine/refine.C src/V_Refine/refine_along_line.C wscript
description:
Experiment with new coarsening and refinement operators. In progress and does not really work.
diff -r cb4e99fc4255 -r b6526d45fde8 src/V_Coarsen/coarsen_2D.C
--- a/src/V_Coarsen/coarsen_2D.C Tue Jun 05 14:14:43 2012 -0700
+++ b/src/V_Coarsen/coarsen_2D.C Tue Jun 05 14:15:49 2012 -0700
@@ -23,6 +23,8 @@
#include "SAMRAI/tbox/Utilities.h"
#include "Constants.h"
+#include "FTensor.hpp"
+
using namespace SAMRAI;
inline void coarsen_point_2D(const pdat::SideIndex &coarse,
@@ -44,9 +46,9 @@ void SAMRAI::geom::V_Coarsen::coarsen_2D
const hier::Box& coarse_box,
const hier::IntVector& ratio) const
{
- const tbox::Dimension& dim(getDim());
+ const tbox::Dimension& Dim(getDim());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, coarse, fine, coarse_box, ratio);
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(Dim, coarse, fine, coarse_box, ratio);
tbox::Pointer<pdat::SideData<double> >
v_fine = fine.getPatchData(src_component);
@@ -65,6 +67,7 @@ void SAMRAI::geom::V_Coarsen::coarsen_2D
const tbox::Pointer<CartesianPatchGeometry> cgeom =
coarse.getPatchGeometry();
+ const double *dx=cgeom->getDx();
/* Numbering of v nodes is
@@ -111,34 +114,128 @@ void SAMRAI::geom::V_Coarsen::coarsen_2D
{
pdat::SideIndex coarse(hier::Index(i,j),0,
pdat::SideIndex::Lower);
- pdat::SideIndex center(coarse*2);
+ pdat::SideIndex fine(coarse*2);
if((i==coarse_box.lower(0)
&& cgeom->getTouchesRegularBoundary(0,0))
|| (i==coarse_box.upper(0)+1
&& cgeom->getTouchesRegularBoundary(0,1)))
{
- (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/2;
+ (*v)(coarse)=((*v_fine)(fine) + (*v_fine)(fine+jp))/2;
}
else
{
+ const int axis=0;
+ FTensor::Tensor1<double,3> offset(0,0,0);
+ offset(axis)=dx[axis]/2;
+ FTensor::Tensor1<double,3> xyz(0,0,0);
+ for(int d=0;d<Dim.getValue();++d)
+ xyz(d)=cgeom->getXLower()[d]
+ + dx[d]*(coarse[d]-coarse_box.lower()[d] + 0.5) - offset(d);
+
coarsen_point_2D(coarse,ip,jp,v,v_fine);
+
+
+ tbox::pout << "coarsen "
+ << coarse << " "
+ << xyz(0) << " "
+ << xyz(1) << " "
+ << dx[0] << " "
+ << fine << " "
+ << (*v)(coarse) << " "
+ << (*v_fine)(fine) << " "
+ << (*v_fine)(fine+ip) << " "
+ << (*v_fine)(fine-ip) << " "
+ << (*v_fine)(fine+jp) << " "
+ << (*v_fine)(fine+jp+ip) << " "
+ << (*v_fine)(fine+jp-ip) << " "
+ << "\n";
}
}
if(directions(1) && i!=coarse_box.upper(0)+1)
{
pdat::SideIndex coarse(hier::Index(i,j),1,
pdat::SideIndex::Lower);
- pdat::SideIndex center(coarse*2);
+ pdat::SideIndex fine(coarse*2);
if((j==coarse_box.lower(1)
&& cgeom->getTouchesRegularBoundary(1,0))
|| (j==coarse_box.upper(1)+1
&& cgeom->getTouchesRegularBoundary(1,1)))
{
- (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/2;
+ (*v)(coarse)=((*v_fine)(fine) + (*v_fine)(fine+ip))/2;
}
else
{
- coarsen_point_2D(coarse,jp,ip,v,v_fine);
+ const int axis=1;
+ FTensor::Tensor1<double,3> offset(0,0,0);
+ offset(axis)=dx[axis]/2;
+ FTensor::Tensor1<double,3> xyz(0,0,0);
+ for(int d=0;d<Dim.getValue();++d)
+ xyz(d)=cgeom->getXLower()[d]
+ + dx[d]*(coarse[d]-coarse_box.lower()[d] + 0.5) - offset(d);
+
+ if(xyz(0)-dx[0]<0.5 && xyz(0)+dx[0]>0.5)
+ {
+ /* Interface between coarse and fine+1 */
+ if((xyz(1)+dx[1]/2>0.6 && xyz(1)<0.6)
+ || (xyz(1)+dx[1]/2>0.4 && xyz(1)<0.4))
+ {
+ // tbox::pout << "coarsen m ";
+ (*v)(coarse)=(((*v_fine)(fine) + (*v_fine)(fine+ip))*2
+ + (*v_fine)(fine-jp) + (*v_fine)(fine+ip-jp))/3;
+ }
+ /* Interface between coarse and fine-1 */
+ else if((xyz(1)-dx[1]/2<0.6 && xyz(1)>0.6)
+ || (xyz(1)-dx[1]/2<0.4 && xyz(1)>0.4))
+ {
+ // tbox::pout << "coarsen p ";
+ (*v)(coarse)=(((*v_fine)(fine) + (*v_fine)(fine+ip))*2
+ + (*v_fine)(fine+jp) + (*v_fine)(fine+ip+jp))/3;
+ }
+ else
+ {
+ // tbox::pout << "coarsen z ";
+ coarsen_point_2D(coarse,jp,ip,v,v_fine);
+ }
+ }
+ else
+ {
+ // tbox::pout << "coarsen ";
+ coarsen_point_2D(coarse,jp,ip,v,v_fine);
+ }
+
+ // /* Interface between coarse and fine-1 */
+ // if(xyz(0)+dx[0]/4>0.5 && xyz(0)<0.5)
+ // {
+ // tbox::pout << "coarsen m ";
+ // (*v)(coarse)=(*v_fine)(fine)/2
+ // + ((*v_fine)(fine-jp) + (*v_fine)(fine+jp))/4;
+ // }
+ // /* Interface between coarse and fine+1 */
+ // else if(xyz(0)-dx[0]/4<0.5 && xyz(0)>0.5)
+ // {
+ // tbox::pout << "coarsen p ";
+ // (*v)(coarse)=(*v_fine)(fine+ip)/2
+ // + ((*v_fine)(fine-jp+ip) + (*v_fine)(fine+jp+ip))/4;
+ // }
+ // else
+ // {
+ // tbox::pout << "coarsen ";
+ // coarsen_point_2D(coarse,jp,ip,v,v_fine);
+ // }
+
+ // tbox::pout << coarse << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << dx[0] << " "
+ // << fine << " "
+ // << (*v)(coarse) << " "
+ // << (*v_fine)(fine) << " "
+ // << (*v_fine)(fine+jp) << " "
+ // << (*v_fine)(fine-jp) << " "
+ // << (*v_fine)(fine+ip) << " "
+ // << (*v_fine)(fine+ip+jp) << " "
+ // << (*v_fine)(fine+ip-jp) << " "
+ // << "\n";
}
}
}
diff -r cb4e99fc4255 -r b6526d45fde8 src/V_Refine.C
--- a/src/V_Refine.C Tue Jun 05 14:14:43 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,158 +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 side-centered double data on
- * a Cartesian mesh.
- *
- ************************************************************************/
-
-#ifndef included_geom_V_Refine_C
-#define included_geom_V_Refine_C
-
-#include "V_Refine.h"
-
-#include <float.h>
-#include <math.h>
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/hier/Index.h"
-#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/pdat/SideVariable.h"
-#include "SAMRAI/tbox/Utilities.h"
-#include "SAMRAI/pdat/CellData.h"
-
-void SAMRAI::geom::V_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::SideOverlap* t_overlap =
- dynamic_cast<const pdat::SideOverlap *>(&fine_overlap);
-
- TBOX_ASSERT(t_overlap != NULL);
-
- for(int axis=0; axis<getDim().getValue(); ++axis)
- {
- const hier::BoxList& boxes = t_overlap->getDestinationBoxList(axis);
- for (hier::BoxList::Iterator b(boxes); b; b++)
- {
- refine(fine,coarse,dst_component,src_component,b(),ratio,axis);
- }
- }
-}
-
-void SAMRAI::geom::V_Refine::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 int &axis) const
-{
- const tbox::Dimension& dimension(getDim());
- const int dim(dimension.getValue());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dimension, fine, coarse, fine_box, ratio);
-
- tbox::Pointer<pdat::SideData<double> >
- v_ptr = coarse.getPatchData(src_component);
- pdat::SideData<double> &v(*v_ptr);
- tbox::Pointer<pdat::SideData<double> >
- v_fine_ptr = fine.getPatchData(dst_component);
- pdat::SideData<double> &v_fine(*v_fine_ptr);
-
-#ifdef DEBUG_CHECK_ASSERTIONS
- TBOX_ASSERT(!v_ptr.isNull());
- TBOX_ASSERT(!v_fine_ptr.isNull());
- TBOX_ASSERT(v.getDepth() == v_fine.getDepth());
- TBOX_ASSERT(v.getDepth() == 1);
-#endif
-
- hier::Box coarse_box=coarse.getBox();
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = coarse.getPatchGeometry();
-
- hier::Index ip(hier::Index::getZeroIndex(dimension)), jp(ip), kp(ip);
- ip[0]=1;
- jp[1]=1;
- if(dim>2)
- kp[2]=1;
- hier::Index pp[]={ip,jp,kp};
-
- for(pdat::CellIterator ci(fine_box); ci; ci++)
- {
- pdat::SideIndex fine(*ci,axis,pdat::SideIndex::Lower);
-
- pdat::SideIndex center(fine);
- center.coarsen(hier::Index::getOneIndex(dimension)*2);
-
- /* This assumes that the levels are always properly nested, so
- that we always have an extra grid space for interpolation.
- So we only have to have a special case for physical
- boundaries, where we do not have an extra grid space. */
-
- double dvx_dy;
-
- if(fine[axis]%2==0)
- {
- /* Maybe this has to be fixed when dvx/dy != 0 on the
- outer boundary because the approximation to the
- derivative is not accurate enough? */
-
- /* Maybe in 3D we should include cross derivatives? */
-
- v_fine(fine)=v(center);
-
- for(int d=(axis+1)%dim;d!=axis;d=(d+1)%dim)
- {
- if(center[d]==coarse_box.lower(d)
- && geom->getTouchesRegularBoundary(d,0))
- {
- dvx_dy=(v(center+pp[d])-v(center))/4;
- }
- else if(center[d]==coarse_box.upper(d)
- && geom->getTouchesRegularBoundary(d,1))
- {
- dvx_dy=(v(center)-v(center-pp[d]))/4;
- }
- else
- {
- dvx_dy=(v(center+pp[d])-v(center-pp[d]))/8;
- }
- v_fine(fine)+=((fine[d]%2==0) ? (-dvx_dy) : dvx_dy);
- }
- }
- else
- {
- v_fine(fine)=(v(center) + v(center+pp[axis]))/2;
-
- for(int d=(axis+1)%dim;d!=axis;d=(d+1)%dim)
- {
- if(center[d]==coarse_box.lower(d)
- && geom->getTouchesRegularBoundary(d,0))
- {
- dvx_dy=(v(center+pp[d])-v(center)
- + v(center+pp[d]+pp[axis])-v(center+pp[axis]))/8;
- }
- else if(center[d]==coarse_box.upper(d)
- && geom->getTouchesRegularBoundary(d,1))
- {
- dvx_dy=(v(center)-v(center-pp[d])
- + v(center+pp[axis])-v(center-pp[d]+pp[axis]))/8;
- }
- else
- {
- dvx_dy=
- (v(center+pp[d])-v(center-pp[d])
- + v(center+pp[d]+pp[axis])-v(center-pp[d]+pp[axis]))/16;
- }
- v_fine(fine)+=((fine[d]%2==0) ? (-dvx_dy) : dvx_dy);
- }
- }
- }
-}
-#endif
diff -r cb4e99fc4255 -r b6526d45fde8 src/V_Refine.h
--- a/src/V_Refine.h Tue Jun 05 14:14:43 2012 -0700
+++ b/src/V_Refine.h Tue Jun 05 14:15:49 2012 -0700
@@ -20,7 +20,9 @@
#include "SAMRAI/hier/Patch.h"
#include "SAMRAI/tbox/Pointer.h"
#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "FTensor.hpp"
#include <string>
namespace SAMRAI {
@@ -133,6 +135,17 @@ public:
const hier::IntVector& ratio,
const int &axis) const;
+ double refine_along_line(pdat::SideData<double> &v,
+ const int &axis,
+ const int &dim,
+ const hier::Index pp[],
+ const pdat::SideIndex &fine,
+ const pdat::SideIndex &coarse,
+ const hier::Box &coarse_box,
+ const CartesianPatchGeometry &coarse_geom,
+ const FTensor::Tensor1<double,3> &xyz,
+ const double *dx) const;
+
private:
std::string d_name_id;
diff -r cb4e99fc4255 -r b6526d45fde8 src/V_Refine/refine.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/V_Refine/refine.C Tue Jun 05 14:15:49 2012 -0700
@@ -0,0 +1,126 @@
+/*************************************************************************
+ *
+ * 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 side-centered double data on
+ * a Cartesian mesh.
+ *
+ ************************************************************************/
+
+#include "V_Refine.h"
+
+#include <float.h>
+#include <math.h>
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/SideVariable.h"
+#include "SAMRAI/tbox/Utilities.h"
+#include "SAMRAI/pdat/CellData.h"
+
+#include "FTensor.hpp"
+
+void SAMRAI::geom::V_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::SideOverlap* t_overlap =
+ dynamic_cast<const pdat::SideOverlap *>(&fine_overlap);
+
+ TBOX_ASSERT(t_overlap != NULL);
+
+ for(int axis=0; axis<getDim().getValue(); ++axis)
+ {
+ const hier::BoxList& boxes = t_overlap->getDestinationBoxList(axis);
+ for (hier::BoxList::Iterator b(boxes); b; b++)
+ {
+ refine(fine,coarse,dst_component,src_component,b(),ratio,axis);
+ }
+ }
+}
+
+void SAMRAI::geom::V_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 int &axis) 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::SideData<double> >
+ v_ptr = coarse_patch.getPatchData(src_component);
+ pdat::SideData<double> &v(*v_ptr);
+ tbox::Pointer<pdat::SideData<double> >
+ v_fine_ptr = fine_patch.getPatchData(dst_component);
+ pdat::SideData<double> &v_fine(*v_fine_ptr);
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!v_ptr.isNull());
+ TBOX_ASSERT(!v_fine_ptr.isNull());
+ TBOX_ASSERT(v.getDepth() == v_fine.getDepth());
+ TBOX_ASSERT(v.getDepth() == 1);
+#endif
+
+ hier::Box coarse_box=coarse_patch.getBox();
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ coarse_geom = coarse_patch.getPatchGeometry();
+
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ fine_geom = fine_patch.getPatchGeometry();
+ const double *dx=fine_geom->getDx();
+
+ const hier::Box &fine_patch_box(fine_patch.getBox());
+
+ hier::Index ip(hier::Index::getZeroIndex(dimension)), jp(ip), kp(ip);
+ ip[0]=1;
+ jp[1]=1;
+ if(dim>2)
+ kp[2]=1;
+ hier::Index pp[]={ip,jp,kp};
+
+ for(pdat::CellIterator ci(fine_box); ci; ci++)
+ {
+ pdat::SideIndex fine(*ci,axis,pdat::SideIndex::Lower);
+
+ pdat::SideIndex coarse(fine);
+ coarse.coarsen(hier::Index::getOneIndex(dimension)*2);
+
+ FTensor::Tensor1<double,3> offset(0,0,0);
+ offset(axis)=dx[axis]/2;
+ FTensor::Tensor1<double,3> xyz(0,0,0);
+ for(int d=0;d<dim;++d)
+ xyz(d)=fine_geom->getXLower()[d]
+ + dx[d]*(fine[d]-fine_patch_box.lower()[d] + 0.5) - offset(d);
+
+ if(fine[axis]%2==0)
+ {
+ v_fine(fine)=
+ refine_along_line(v,axis,dim,pp,fine,coarse,coarse_box,
+ *coarse_geom,xyz,dx);
+ }
+ else
+ {
+ FTensor::Tensor1<double,3> xyz_low, xyz_high;
+ FTensor::Index<'a',3> a;
+
+ xyz_low(a)=xyz(a) - 2*offset(a);
+ xyz_high(a)=xyz(a) + 2*offset(a);
+
+ v_fine(fine)=
+ (refine_along_line(v,axis,dim,pp,fine,coarse,coarse_box,
+ *coarse_geom,xyz_low,dx)
+ + refine_along_line(v,axis,dim,pp,fine,coarse+pp[axis],
+ coarse_box,*coarse_geom,xyz_high,dx))/2;
+ }
+ }
+}
diff -r cb4e99fc4255 -r b6526d45fde8 src/V_Refine/refine_along_line.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/V_Refine/refine_along_line.C Tue Jun 05 14:15:49 2012 -0700
@@ -0,0 +1,280 @@
+#include "V_Refine.h"
+
+/* This assumes that the levels are always properly nested, so that we
+ always have an extra grid space for interpolation. So we only have
+ to have a special case for physical boundaries, where we do not
+ have an extra grid space. */
+
+/* Maybe this has to be fixed when dvx/dy != 0 on the outer boundary
+ because the approximation to the derivative is not accurate
+ enough? */
+
+/* Maybe in 3D we should include cross derivatives? */
+
+double SAMRAI::geom::V_Refine::refine_along_line
+(pdat::SideData<double> &v,
+ const int &axis,
+ const int &dim,
+ const hier::Index pp[],
+ const pdat::SideIndex &fine,
+ const pdat::SideIndex &coarse,
+ const hier::Box &coarse_box,
+ const CartesianPatchGeometry &coarse_geom,
+ const FTensor::Tensor1<double,3> &xyz,
+ const double *dx) const
+{
+ double result=v(coarse);
+
+ for(int d=(axis+1)%dim;d!=axis;d=(d+1)%dim)
+ {
+ const int sgn(fine[d]%2==0 ? -1 : 1);
+
+ double dvx_dy;
+ if(coarse[d]==coarse_box.lower(d)
+ && coarse_geom.getTouchesRegularBoundary(d,0))
+ {
+ dvx_dy=sgn*(v(coarse+pp[d])-v(coarse))/4;
+ }
+ else if(coarse[d]==coarse_box.upper(d)
+ && coarse_geom.getTouchesRegularBoundary(d,1))
+ {
+ dvx_dy=sgn*(v(coarse)-v(coarse-pp[d]))/4;
+ }
+ else
+ {
+ dvx_dy=sgn*(v(coarse+pp[d])-v(coarse-pp[d]))/8;
+
+ if(axis==0 && xyz(0)-dx[0]<0.5+1e-6 && xyz(0)+dx[0]>0.5+1e-6)
+ {
+ const double y_min(0.5-sgn*0.1), y_max(0.5+sgn*0.1);
+ /* Top tip */
+ /* Singularity in coarse+1 */
+ if(sgn*(xyz(1)+sgn*dx[1]*0.5-y_max)<0 && sgn*(xyz(1)+sgn*dx[1]*2.5-y_max)>0)
+ {
+ dvx_dy=(v(coarse)-v(coarse-pp[d]*sgn))/4;
+ // tbox::pout << "refine p cp1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+
+ }
+ /* Singularity in coarse */
+ else if(sgn*(xyz(1)-sgn*dx[1]*1.5-y_max)<0 && sgn*(xyz(1)+sgn*dx[1]*0.5-y_max)>0)
+ {
+ dvx_dy=0;
+ // tbox::pout << "refine p c "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ }
+ /* Singularity in coarse-1 */
+ else if(sgn*(xyz(1)-sgn*dx[1]*3.5-y_max)<0 && sgn*(xyz(1)-sgn*dx[1]*1.5-y_max)>0)
+ {
+ dvx_dy=(v(coarse+pp[d]*sgn)-v(coarse))/4;
+ // tbox::pout << "refine p cm1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ }
+
+ /* Bottom tip */
+ /* Singularity in coarse+1 */
+ else if(sgn*(xyz(1)+sgn*dx[1]*0.5-y_min)<0 && sgn*(xyz(1)+sgn*dx[1]*2.5-y_min)>0)
+ {
+ dvx_dy=(v(coarse)-v(coarse-pp[d]*sgn))/4;
+ // tbox::pout << "refine m cp1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ }
+ /* Singularity in coarse */
+ else if(sgn*(xyz(1)-sgn*dx[1]*1.5-y_min)<0 && sgn*(xyz(1)+sgn*dx[1]*0.5-y_min)>0)
+ {
+ dvx_dy=0;
+ // tbox::pout << "refine m c "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ }
+ /* Singularity in coarse-1 */
+ else if(sgn*(xyz(1)-sgn*dx[1]*3.5-y_min)<0 && sgn*(xyz(1)-sgn*dx[1]*1.5-y_min)>0)
+ {
+ dvx_dy=(v(coarse+pp[d]*sgn)-v(coarse))/4;
+ // tbox::pout << "refine m cm1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ }
+
+ // /* Bottom tip */
+ // /* Singularity in coarse+1 */
+ // else if(sgn*(xyz(1)-sgn*dx[1]*1.5-y_min)<0 && sgn*(xyz(1)-sgn*dx[1]*3.5-y_min)>0)
+ // {
+ // dvx_dy=(v(coarse)-v(coarse-pp[d]*sgn))/4;
+ // tbox::pout << "refine m cp1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ // }
+ // /* Jump in coarse */
+ // else if(sgn*(xyz(1)-sgn*dx[1]*1.5-y_min)>0 && sgn*(xyz(1)+sgn*dx[1]*0.5-y_min)<0)
+ // {
+ // dvx_dy=(v(coarse+pp[d]*sgn)-v(coarse))/4;
+ // tbox::pout << "refine m c "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+ // }
+
+ }
+
+ // if(coarse[1]==6)
+ // {
+ // tbox::pout << "refine "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + dvx_dy << " "
+ // << "\n";
+
+ // }
+
+
+
+ if(axis==1 && xyz(1)>=0.4 && xyz(1)<=0.6)
+ {
+ /* Interface between fine and coarse+1, do a one-sided
+ interpolation. */
+ if(sgn*(xyz(0)-0.5)<0 && sgn*(xyz(0)+sgn*1.5*dx[0]-0.5)>0)
+ {
+ dvx_dy=(v(coarse)-v(coarse-pp[d]*sgn))/4;
+
+ // tbox::pout << "refine fcp1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + sgn*dvx_dy << " "
+ // << "\n";
+ }
+ /* Interface between fine and coarse, use derivative
+ from the other side. */
+ else if(sgn*(xyz(0)-0.5)>0 && sgn*(xyz(0)-sgn*dx[0]/2-0.5)<0)
+ {
+ dvx_dy=(v(coarse+pp[d]*sgn)-v(coarse)
+ - (v(coarse)-v(coarse-pp[d]*sgn))*0.75);
+
+ // tbox::pout << "refine fc "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + sgn*dvx_dy << " "
+ // << "\n";
+ }
+ /* Interface between coarse and coarse-1, do a one-sided
+ interpolation. */
+ else if(sgn*(xyz(0)-0.5)>0 && sgn*(xyz(0)-sgn*2.5*dx[0]-0.5)<0)
+ {
+ dvx_dy=(v(coarse+pp[d]*sgn)-v(coarse))/4;
+
+ // tbox::pout << "refine ccm1 "
+ // << 1/dx[0] << " "
+ // << xyz(0) << " "
+ // << xyz(1) << " "
+ // << fine << " "
+ // << coarse << " "
+ // << sgn << " "
+ // << v(coarse+pp[d]*sgn) << " "
+ // << v(coarse) << " "
+ // << v(coarse-pp[d]*sgn) << " "
+ // << result + sgn*dvx_dy << " "
+ // << "\n";
+ }
+ }
+ }
+ result+=dvx_dy;
+ }
+ return result;
+}
+
diff -r cb4e99fc4255 -r b6526d45fde8 wscript
--- a/wscript Tue Jun 05 14:14:43 2012 -0700
+++ b/wscript Tue Jun 05 14:15:49 2012 -0700
@@ -7,7 +7,6 @@ def configure(conf):
def configure(conf):
conf.load('compiler_cxx')
def build(bld):
- # bld.recurse('src/libtcod-1.5.0')
bld.program(
features = ['cxx','cprogram'],
source = ['src/main.C',
@@ -20,7 +19,8 @@ def build(bld):
'src/FACStokes/setupPlotter.C',
'src/FACStokes/solveStokes.C',
'src/P_Refine.C',
- 'src/V_Refine.C',
+ 'src/V_Refine/refine.C',
+ 'src/V_Refine/refine_along_line.C',
'src/Resid_Coarsen.C',
'src/V_Coarsen/coarsen_2D.C',
'src/V_Coarsen/coarsen_3D.C',
@@ -51,6 +51,7 @@ def build(bld):
'src/StokesFACOps/smooth_Tackley_3D.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',
@@ -73,26 +74,44 @@ def build(bld):
'src/StokesHypreSolver.C'],
target = 'gamr',
- # cxxflags = ['-std=c++0x','-g'],
- cxxflags = ['-O3', '-Wall', '-Wextra', '-Wconversion',
- '-DTESTING=0'],
+ # cxxflags = ['-O3', '-Wall', '-Wextra', '-Wconversion',
+ # '-DTESTING=0', '-Drestrict='],
+ cxxflags = ['-g', '-Wall', '-Wextra', '-Wconversion',
+ '-Drestrict='],
lib = ['SAMRAI_appu', 'SAMRAI_algs', 'SAMRAI_solv',
'SAMRAI_geom', 'SAMRAI_mesh', 'SAMRAI_math',
'SAMRAI_pdat', 'SAMRAI_xfer', 'SAMRAI_hier',
'SAMRAI_tbox', 'petsc', 'hdf5',
+ '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',
'dl', 'm', 'gfortranbegin', 'gfortran', 'm',
'gfortranbegin', 'gfortran', 'm', 'gfortranbegin',
'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',
+ libpath = ['/home/boo/cig/SAMRAI-hg-build/lib',
+ '/usr/lib/petsc/linux-gnu-c-opt/lib'],
+ # 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','src/FTensor','../../SAMRAI-hg-build/include',
+ '../../SAMRAI-hg/source',
+ '/usr/lib/petsc/include',
'/usr/lib/petsc/bmake/linux-gnu-c-opt',
- '../../petsc/petsc-3.2-p7/src/vec',
+ '/usr/lib/petsc/src/vec',
'/usr/include']
+ # includes = ['src','../../samrai/include',
+ # '../../samrai/source',
+ # '../../petsc/petsc-3.2-p7/real8-with-debug/include',
+ # '/usr/lib/petsc/bmake/linux-gnu-c-opt',
+ # '../../petsc/petsc-3.2-p7/src/vec',
+ # '/usr/include']
)
More information about the CIG-COMMITS
mailing list