testdivurm.cpp

强大的并行科学计算软件包

// -*- C++ -*-
// ACL:license
// ----------------------------------------------------------------------
// This software and ancillary information (herein called "SOFTWARE")
// called POOMA (Parallel Object-Oriented Methods and Applications) is
// made available under the terms described here.  The SOFTWARE has been
// approved for release with associated LA-CC Number LA-CC-98-65.
// 
// Unless otherwise indicated, this SOFTWARE has been authored by an
// employee or employees of the University of California, operator of the
// Los Alamos National Laboratory under Contract No.  W-7405-ENG-36 with
// the U.S. Department of Energy.  The U.S. Government has rights to use,
// reproduce, and distribute this SOFTWARE, and to allow others to do so.
// The public may copy and use this SOFTWARE, FOR NONCOMMERCIAL USE ONLY,
// without charge, provided that this Notice and any statement of
// authorship are reproduced on all copies.  Neither the Government nor
// the University makes any warranty, express or implied, or assumes any
// liability or responsibility for the use of this SOFTWARE.
// 
// If SOFTWARE is modified to produce derivative works, such modified
// SOFTWARE should be clearly marked, so as not to confuse it with the
// version available from LANL.
// 
// For more information about POOMA, send e-mail to [e-mail],
// or visit the POOMA web page at http://www.acl.lanl.gov/pooma/.
// ----------------------------------------------------------------------
// ACL:license

//-----------------------------------------------------------------------------
// TestDivURM: Test Divergence function with a Uniform Rectilinear Mesh.
//-----------------------------------------------------------------------------

#include "Pooma/Fields.h"
#include "Utilities/Tester.h"

// Forward declarations:

const double epsilon = 1.0e-6;

// Vert-->Cell

template<int Dim>
void doTestVC(Field<DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >, 
  Vector<Dim>, MultiPatch<GridTag, Brick> > &fv, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Div<Cell, 
    DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >, Vector<Dim> > >,
    ConstField<DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >,
    Vector<Dim>, MultiPatch<GridTag, Brick> > >::Type_t d = div<Cell>(fv);
    
  Interval<Dim> outCells = d.physicalDomain();
  int nCellsOut = outCells.size();
  double error = abs(sum(d)/nCellsOut - 1.0*Dim);
  if (error > epsilon) {
    tester.out() << "div<Cell>(" << Dim << "D Vert position field) != const="
	      << 1.0*Dim << " ; sum = " 
	      << error
	      << std::endl;
    tester.check(false);
    }
}

// Cell-->Vert

template<int Dim>
void doTestCV(Field<DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >, 
  Vector<Dim>, MultiPatch<GridTag, Brick> > &fc, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Div<Vert, 
    DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >, Vector<Dim> > >,
    ConstField<DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >,
    Vector<Dim>, MultiPatch<GridTag, Brick> > >::Type_t d = div<Vert>(fc);
    
  Interval<Dim> outVerts = d.physicalDomain();
  int nVertsOut = outVerts.size();
  double error = abs(sum(d)/nVertsOut - 1.0*Dim);
  if (error > epsilon) {
    tester.out() << "div<Vert>(" << Dim << "D Cell position field) != const="
	      << 1.0*Dim << " ; sum = " 
	      << error
	      << std::endl;
    tester.check(false);
    }
}

// Vert-->Vert

template<int Dim>
void doTestVV(Field<DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >, 
  Vector<Dim>, MultiPatch<GridTag, Brick> > &fv, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Div<Vert, 
    DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >, Vector<Dim> > >,
    ConstField<DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >,
    Vector<Dim>, MultiPatch<GridTag, Brick> > >::Type_t d = div<Vert>(fv);
    
  Interval<Dim> outVerts = d.physicalDomain();
  int nVertsOut = outVerts.size();
  double error = abs(sum(d)/nVertsOut - 1.0*Dim);
  if (error > epsilon) {
    tester.out() << "div<Vert>(" << Dim << "D Vert position field) != const="
	      << 1.0*Dim << " ; sum = " 
	      << error
	      << std::endl;
    tester.check(false);
    }
}

// Cell-->Cell

template<int Dim>
void doTestCC(Field<DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >, 
  Vector<Dim>, MultiPatch<GridTag, Brick> > &fc, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Div<Cell, 
    DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >, Vector<Dim> > >,
    ConstField<DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >,
    Vector<Dim>, MultiPatch<GridTag, Brick> > >::Type_t d = div<Cell>(fc);
    
  Interval<Dim> outCells = d.physicalDomain();
  int nCellsOut = outCells.size();
  double error = abs(sum(d)/nCellsOut - 1.0*Dim);
  if (error > epsilon) {
    tester.out() << "div<Cell>(" << Dim << "D Cell position field) != const="
	      << 1.0*Dim << " ; sum = " 
	      << error
	      << std::endl;
    tester.check(false);
    }
}

// Tensor -> Vector, Vert-->Cell

template<int Dim>
void doTestTVVC(Field<DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >, 
  Tensor<Dim>, MultiPatch<GridTag, Brick> > &ftv, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Div<Cell, 
    DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >, Tensor<Dim> > >,
    ConstField<DiscreteGeometry<Vert, UniformRectilinearMesh<Dim> >,
    Tensor<Dim>, MultiPatch<GridTag, Brick> > >::Type_t d = div<Cell>(ftv);
    
  Interval<Dim> outCells = d.physicalDomain();
  int nCellsOut = outCells.size();
  double summer = sum(dot(d, d)) / nCellsOut;
  Vector<Dim,double> dees(1.0*Dim);
  double comparisonVectorMag = dot(dees, dees);
  if (abs(abs(summer) - comparisonVectorMag) > epsilon) {
    tester.out() << "div<Cell>(" << Dim 
	      << "D Vert position Tensor field) != const="
	      << comparisonVectorMag << " ; summer = " << summer 
	      << std::endl;
    tester.check(false);
  }
}

// Tensor -> Vector, Cell-->Vert

template<int Dim>
void doTestTVCV(Field<DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >, 
  Tensor<Dim>, MultiPatch<GridTag, Brick> > &ftc, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Div<Vert, 
    DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >, Tensor<Dim> > >,
    ConstField<DiscreteGeometry<Cell, UniformRectilinearMesh<Dim> >,
    Tensor<Dim>, MultiPatch<GridTag, Brick> > >::Type_t d = div<Vert>(ftc);
    
  Interval<Dim> outVerts = d.physicalDomain();
  int nVertsOut = outVerts.size();
  double summer = sum(dot(d, d)) / nVertsOut;
  Vector<Dim,double> dees(1.0*Dim);
  double comparisonVectorMag = dot(dees, dees);
  if (abs(abs(summer) - comparisonVectorMag) > epsilon) {
    tester.out() << "div<Vert>(" << Dim 
	      << "D Vert position Tensor field) != const="
	      << comparisonVectorMag << " ; summer = " << summer 
	      << std::endl;
    tester.check(false);
  }
}

template<int Dim>
void divTest(const int *nVerts, Pooma::Tester &tester)
{
  int i, d, d2, d3;
  double summer = 0.0;
  int nCellsOut, nVertsOut;

  // Vertex and cell domains:
  Interval<Dim> vertDomain, cellDomain;
  for (d = 0; d < Dim; d++) {
    vertDomain[d] = Interval<1>(nVerts[d]);
    cellDomain[d] = Interval<1>(nVerts[d] - 1);
  }

  // Create the mesh:
  typedef UniformRectilinearMesh<Dim> Mesh_t;
  Mesh_t::PointType_t origin(0.0);
  Vector<Dim,double> spacings;
  for (d = 0; d < Dim; d++) {
    origin(d) = d;
    spacings(d) = d + 1;
  }
  Mesh_t mesh(vertDomain, origin, spacings);
  
  // Create the geometries:
  DiscreteGeometry<Cell, Mesh_t> geomc(mesh, GuardLayers<Dim>(1));
  DiscreteGeometry<Vert, Mesh_t> geomv(mesh, GuardLayers<Dim>(1));
  
  // Make the fields:
  Loc<Dim> blocks(2);
  GridLayout<Dim> layoutc(cellDomain, blocks, 
                          GuardLayers<Dim>(0), GuardLayers<Dim>(1));
  GridLayout<Dim> layoutv(vertDomain, blocks, 
                          GuardLayers<Dim>(0), GuardLayers<Dim>(1));
  typedef MultiPatch<GridTag, Brick> MP_t;
  Field<DiscreteGeometry<Cell, Mesh_t >, Vector<Dim>, MP_t> fc(geomc, layoutc);
  Field<DiscreteGeometry<Vert, Mesh_t >, Vector<Dim>, MP_t> fv(geomv, layoutv);
  Field<DiscreteGeometry<Vert, Mesh_t >, Tensor<Dim>, MP_t> ftv(geomv, layoutv);
  Field<DiscreteGeometry<Cell, Mesh_t >, Tensor<Dim>, MP_t> ftc(geomc, layoutc);

  // Set the Vector Field values equal to the positions
  Interval<Dim> tdv = fv.totalDomain();
  fv(tdv) = fv.x(tdv);
  Interval<Dim> tdc = fc.totalDomain();
  fc(tdc) = fc.x(tdc);

  // Assign positive-sloping, linear-ramp values into the Tensor Field:
  ftv(tdv) = Tensor<Dim>(0.0);
  for (d = 0; d < Dim; d++ ) {
    for (d2 = 0; d2 < Dim; d2++) {
      for (d3 = 0; d3 < Dim; d3++) {
	ftv.comp(d,d2)(tdv) += fv.comp(d3)(tdv);
      }
    }
  }
  ftc(tdc) = Tensor<Dim>(0.0);
  for (d = 0; d < Dim; d++ ) {
    for (d2 = 0; d2 < Dim; d2++) {
      for (d3 = 0; d3 < Dim; d3++) {
	ftc.comp(d,d2)(tdc) += fc.comp(d3)(tdc);
      }
    }
  }
      
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // scalar->Vector:
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

  // Vert-->Cell
  doTestVC(fv, tester);

  // Cell-->Vert
  doTestCV(fc, tester);

  // Vert-->Vert
  doTestVV(fv, tester);

  // Cell-->Cell
  doTestCC(fc, tester);

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // Tensor->Vector:
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

  // Vert-->Cell
  doTestTVVC(ftv, tester);

  // Cell-->Vert
  doTestTVCV(ftc, tester);
}
   
int main(int argc, char *argv[])
{
  Pooma::initialize(argc,argv);
  Pooma::Tester tester(argc, argv);

  // Create the sizes of the physical (vertex) domains (1D, 2D, and 3D):
  const int nVerts[3] = {4, 4, 4};

  // Call the test subroutine, specifying template instances for {1,2,3}D:

  // If compiling with CW, only do a subset of the tests.
  
  divTest<1>(nVerts, tester); // 1D

#if !defined(__MWERKS__)

  divTest<2>(nVerts, tester); // 2D
  divTest<3>(nVerts, tester); // 3D
  
#endif

  int ret = tester.results("TestDivURM");
  Pooma::finalize();
  return ret;
}

// ACL:rcsinfo
// ----------------------------------------------------------------------
// $RCSfile: TestDivURM.cpp,v $   $Author: swhaney $
// $Revision: 1.8 $   $Date: 1999/08/20 18:57:31 $
// ----------------------------------------------------------------------
// ACL:rcsinfo