testgradrm.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

//-----------------------------------------------------------------------------
// TestGradRM: Test Gradient function with a 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, RectilinearMesh<Dim> >, 
  double, MultiPatch<GridTag, Brick> > &fsv, Pooma::Tester &tester)
{
  typename View1<FieldStencil<Grad<Cell, 
    DiscreteGeometry<Vert, RectilinearMesh<Dim> >, double> >,
    ConstField<DiscreteGeometry<Vert, RectilinearMesh<Dim> >,
    double, MultiPatch<GridTag, Brick> > >::Type_t g = grad<Cell>(fsv);
    
  Interval<Dim> outCells = g.physicalDomain();
  int nCellsOut = outCells.size();
  double summer = sum(dot(g, g)) / nCellsOut;
  double comparisonVectorMag = Dim;
  if (abs(abs(summer) - comparisonVectorMag) > epsilon) {
    tester.out() << "grad<Cell>(" << Dim 
	      << "D Vert position scalar field) != const="
	      << comparisonVectorMag << ") ; summer = " << summer << std::endl;
    tester.check(false);
  }
}

// Cell-->Vert

template<int Dim>
void doTestCV(Field<DiscreteGeometry<Cell, RectilinearMesh<Dim> >, 
  double, MultiPatch<GridTag, Brick> > &fsc, Pooma::Tester &tester)
{
  Interval<Dim> outVerts = grad<Vert>(fsc).physicalDomain();
  int nVertsOut = outVerts.size();
  double summer = 
    sum(dot(grad<Vert>(fsc)(outVerts), grad<Vert>(fsc)(outVerts))) / nVertsOut;
  Vector<Dim,double> ones(1.0);
  double comparisonVectorMag = dot(ones, ones);
  if (abs(abs(summer) - comparisonVectorMag) > epsilon) {
    tester.out() << "grad<Vert>(" << Dim 
	      << "D Cell position scalar field) != const="
	      << comparisonVectorMag << ") ; summer = " << summer << std::endl;
    tester.check(false);
  }
}

// Vert-->Vert

template<int Dim>
void doTestVV(Field<DiscreteGeometry<Vert, RectilinearMesh<Dim> >, 
  double, MultiPatch<GridTag, Brick> > &fsv, Pooma::Tester &tester)
{
  Interval<Dim> outVerts = grad<Vert>(fsv).physicalDomain();
  int nVertsOut = outVerts.size();
  double summer = 
    sum(dot(grad<Vert>(fsv)(outVerts), grad<Vert>(fsv)(outVerts))) / nVertsOut;
  Vector<Dim,double> ones(1.0);
  double comparisonVectorMag = dot(ones, ones);
  if (abs(abs(summer) - comparisonVectorMag) > epsilon) {
    tester.out() << "grad<Vert>(" << Dim 
	      << "D Vert position scalar field) != const="
	      << comparisonVectorMag << ") ; summer = " << summer << std::endl;
    tester.check(false);
  }
}

// Cell-->Cell

template<int Dim>
void doTestCC(Field<DiscreteGeometry<Cell, RectilinearMesh<Dim> >, 
  double, MultiPatch<GridTag, Brick> > &fsc, Pooma::Tester &tester)
{
  Interval<Dim> outCells = grad<Cell>(fsc).physicalDomain();
  int nCellsOut = outCells.size();
  double summer = 
    sum(dot(grad<Cell>(fsc)(outCells), grad<Cell>(fsc)(outCells))) / nCellsOut;
  Vector<Dim,double> ones(1.0);
  double comparisonVectorMag = dot(ones, ones);
  if (abs(abs(summer) - comparisonVectorMag) > epsilon) {
    tester.out() << "grad<Cell>(" << Dim 
	      << "D Cell position scalar field) != const="
	      << comparisonVectorMag << ") ; summer = " << summer << std::endl;
    tester.check(false);
  }
}

// Vector -> Tensor, Vert-->Cell
// Since this is the gradient of the position vector (x*x_hat), the result
// should be the identity tensor (NRL Plasma Formulary Vector Identities
// section):

template<int Dim>
void doTestVTVC(Field<DiscreteGeometry<Vert, RectilinearMesh<Dim> >, 
  Vector<Dim>, MultiPatch<GridTag, Brick> > &fv, Pooma::Tester &tester)
{
  Tensor<Dim,double> identityt(0.0);
  int d, d2;
  for (d = 0; d < Dim; d++) { identityt(d,d) = 1.0; }
  Interval<Dim> outCells = grad<Cell>(fv).physicalDomain();
  int nCellsOut = outCells.size();
  Tensor<Dim,double> summert(0.0);
  summert = 
    sum( sqrt( (grad<Cell>(fv)(outCells) - identityt) * 
               (grad<Cell>(fv)(outCells) - identityt) ) ) / nCellsOut;
  double comparisonTensorMag = 0.0;
  for (d = 0; d < Dim; d++) {
    for (d2 = 0; d2 < Dim; d2++) {
      comparisonTensorMag += summert(d,d2);
    }
  }
  if (comparisonTensorMag > epsilon) {
    tester.out() << "grad<Cell>(" << Dim
	      << "D Vert position vector field) != const="
	      << identityt << ") ; summert = " << summert << std::endl;
    tester.check(false);
  }
}

// Vector -> Tensor, Cell-->Vert
// Since this is the gradient of the position vector (x*x_hat), the result
// should be the identity tensor (NRL Plasma Formulary Vector Identities
// section):

template<int Dim>
void doTestVTCV(Field<DiscreteGeometry<Cell, RectilinearMesh<Dim> >, 
  Vector<Dim>, MultiPatch<GridTag, Brick> > &fc, Pooma::Tester &tester)
{
  Tensor<Dim,double> identityt(0.0);
  int d, d2;
  for (d = 0; d < Dim; d++) { identityt(d,d) = 1.0; }
  Interval<Dim> outVerts = grad<Vert>(fc).physicalDomain();
  int nVertsOut = outVerts.size();
  Tensor<Dim,double> summert(0.0);
  summert = 
    sum( sqrt( (grad<Vert>(fc)(outVerts) - identityt) * 
               (grad<Vert>(fc)(outVerts) - identityt) ) ) / nVertsOut;
  double comparisonTensorMag = 0.0;
  for (d = 0; d < Dim; d++) {
    for (d2 = 0; d2 < Dim; d2++) {
      comparisonTensorMag += summert(d,d2);
    }
  }
  if (comparisonTensorMag > epsilon) {
    tester.out() << "grad<Vert>(" << Dim
	      << "D Cell position vector field) != const="
	      << identityt << ") ; summert = " << summert << std::endl;
    tester.check(false);
  }
}

template<int Dim>
void gradTest(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 RectilinearMesh<Dim> Mesh_t;
  Mesh_t::PointType_t origin(0.0);
  Vector<Dim, Array<1, Mesh_t::AxisType_t> > spacings;
  for (d = 0; d < Dim; d++) {
    spacings(d).initialize(cellDomain[d]);
    for (i = 0; i < cellDomain[d].size(); i++) {
      spacings(d)(i) = origin(d) + (d + 2)*(i + 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 >, double, MP_t> fsc(geomc, layoutc);
  Field<DiscreteGeometry<Vert, Mesh_t >, double, MP_t> fsv(geomv, layoutv);
  Field<DiscreteGeometry<Vert, Mesh_t >, Vector<Dim>, MP_t> fv(geomv, layoutv);
  Field<DiscreteGeometry<Cell, Mesh_t >, Vector<Dim>, MP_t> fc(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 scalar Fields:
  // Fields:
  fsv(tdv) = 0.0;
  for (d = 0; d < Dim; d++ ) {
    fsv(tdv) += fv.comp(d)(tdv);
  }
  fsc(tdc) = 0.0;
  for (d = 0; d < Dim; d++ ) {
    fsc(tdc) += fc.comp(d)(tdc);
  }
      
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  // scalar->Vector:
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

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

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

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

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

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

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

  // Cell-->Vert
  doTestVTCV(fc, 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.
  
  gradTest<1>(nVerts, tester); // 1D

#if !defined(__MWERKS__)

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

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

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