lux_test2.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
//-----------------------------------------------------------------------------
// Lux test 2: Display a 2D and a 3D particle set
//-----------------------------------------------------------------------------

#include "Pooma/Pooma.h"
#include "Pooma/Lux.h"
#include "Pooma/Arrays.h"
#include "Pooma/Fields.h"
#include "Pooma/Particles.h"
#include "Pooma/Domains.h"
#include "Utilities/Tester.h"


// Traits class for Particles object
template <class EngineTag, class Centering, class MeshType, class FL,
          class InterpolatorTag>
struct PTraits
{
  // The type of engine to use in the attributes
  typedef EngineTag AttributeEngineTag_t;

  // The type of particle layout to use
  typedef SpatialLayout<DiscreteGeometry<Centering,MeshType>,FL> 
    ParticleLayout_t;

  // The type of interpolator to use
  typedef InterpolatorTag InterpolatorTag_t;
};

// Particles subclass with position and velocity
template <class PT>
class ChargedParticles : public Particles<PT>
{
public:
  // Typedefs
  typedef Particles<PT>                          Base_t;
  typedef typename Base_t::AttributeEngineTag_t  AttributeEngineTag_t;
  typedef typename Base_t::ParticleLayout_t      ParticleLayout_t;
  typedef typename ParticleLayout_t::AxisType_t  AxisType_t;
  typedef typename ParticleLayout_t::PointType_t PointType_t;
  typedef typename PT::InterpolatorTag_t         InterpolatorTag_t;

  // Dimensionality
  static const int dimensions = ParticleLayout_t::dimensions;

  // Constructor: set up layouts, register attributes
  ChargedParticles(const ParticleLayout_t &pl)
  : Particles<PT>(pl)
  {
    addAttribute(R);
    addAttribute(V);
    addAttribute(E);
    addAttribute(qm);
  }

  // Position and velocity attributes (as public members)
  DynamicArray<PointType_t,AttributeEngineTag_t> R;
  DynamicArray<PointType_t,AttributeEngineTag_t> V;
  DynamicArray<PointType_t,AttributeEngineTag_t> E;
  DynamicArray<double,     AttributeEngineTag_t> qm;
};


// Dimensionality of this problem
static const int PDim = 2;

// Engine tag type for attributes
typedef MultiPatch<GridTag,Brick> AttrEngineTag_t;

// Mesh type
typedef UniformRectilinearMesh<PDim,Cartesian<PDim>,double> Mesh_t;

// Centering of Field elements on mesh
typedef Cell Centering_t;

// Geometry type for Fields
typedef DiscreteGeometry<Centering_t,Mesh_t> Geometry_t;

// Field types
typedef Field< Geometry_t, double,
               MultiPatch<UniformTag,Brick> > DField_t;
typedef Field< Geometry_t, Vector<PDim,double>,
               MultiPatch<UniformTag,Brick> > VecField_t;

// Field layout type, derived from Engine type
typedef DField_t::Engine_t Engine_t;
typedef Engine_t::Layout_t FLayout_t;

// Type of interpolator
typedef NGP InterpolatorTag_t;

// Particle traits class
typedef PTraits<AttrEngineTag_t,Centering_t,Mesh_t,FLayout_t,
                InterpolatorTag_t> PTraits_t;

// Type of particle layout
typedef PTraits_t::ParticleLayout_t PLayout_t;

// Type of actual particles
typedef ChargedParticles<PTraits_t> Particles_t;

// Grid sizes
const int nx = 200, ny = 200;

// Number of particles in simulation
const int NumPart = 400;

// Number of timesteps in simulation
const int NumSteps = 20;

// The value of pi
const double pi = acos(-1.0);

// Maximum value for particle q/m ratio
const double qmmax = 1.0;

// Timestep
const double dt = 1.0;


int main(int argc, char *argv[])
{
  // Initialize POOMA and output stream, using Tester class

  Pooma::initialize(argc, argv);
  Pooma::Tester tester(argc, argv);
  tester.out() << argv[0] << ": Lux Particles PIC2d display test" << std::endl;
  tester.out() << "---------------------------------------------" << std::endl;

  tester.out() << "Initializing particles ..." << std::endl;

  // Create mesh and geometry objects for cell-centered fields.
  Interval<PDim> meshDomain(nx+1,ny+1);
  Mesh_t mesh(meshDomain);
  Geometry_t geometry(mesh);

  // Create a second geometry object that includes a guard layer.
  GuardLayers<PDim> gl(1);
  Geometry_t geometryGL(mesh,gl);

  // Create field layout objects for our electrostatic potential
  // and our electric field.  Decomposition is 4 x 4.
  Loc<PDim> blocks(4,4);
  FLayout_t flayout(geometry.physicalDomain(),blocks);
  FLayout_t flayoutGL(geometryGL.physicalDomain(),blocks,gl);

  // Create and initialize electrostatic potential and electric field.
  DField_t phi(geometryGL,flayoutGL);
  VecField_t EFD(geometry,flayout);

  // potential phi = phi0 * sin(2*pi*x/Lx) * cos(4*pi*y/Ly)
  // Note that phi is a periodic Field
  // Electric field EFD = -grad(phi);
  phi.addBoundaryConditions(AllPeriodicFaceBC());
  double phi0 = 0.01 * static_cast<double>(nx);
  phi = phi0 * sin(2.0*pi*phi.x().comp(0)/nx)
             * cos(4.0*pi*phi.x().comp(1)/ny);
  EFD = -grad<Centering_t>(phi);

  // Create a particle layout object for our use
  PLayout_t layout(geometry,flayout);

  // Create a Particles object and set periodic boundary conditions
  Particles_t P(layout);
  Particles_t::PointType_t lower(0.0,0.0), upper(nx,ny);
  PeriodicBC<Particles_t::PointType_t> bc(lower,upper);
  P.addBoundaryCondition(P.R,bc);

  // Create an equal number of particles on each processor
  // and recompute global domain.
  P.globalCreate(NumPart);

  // Random initialization for particle positions in nx by ny domain
  // Zero initialization for particle velocities
  // Random intialization for charge-to-mass ratio from -qmmax to qmmax
  P.V = Particles_t::PointType_t(0.0);
  srand(12345U);
  Particles_t::PointType_t initPos;
  for (int i = 0; i < NumPart; ++i)
  {
    initPos(0) = nx * rand() /
	         static_cast<Particles_t::AxisType_t>(RAND_MAX);
    initPos(1) = ny * rand() /
	         static_cast<Particles_t::AxisType_t>(RAND_MAX);
    P.R(i) = initPos;
    P.qm(i) = (2.0 * rand() / static_cast<double>(RAND_MAX) - 1.0) *
              qmmax;
  }

  // Redistribute particle data based on spatial layout
  P.swap(P.R);

  tester.out() << "PIC2d setup complete." << std::endl;
  tester.out() << "---------------------" << std::endl;

  // Create a Lux connection

  tester.out() << "Creating LuxConnection object ..." << std::endl;
  Connection<Lux> lux("test2");
  tester.out() << "Finished creating LuxConnection object." << std::endl;

  // Add attributes in to the display

  tester.out() << "Connecting qm for display ..." << std::endl;
  lux.connect("P-qm", P.R, P.qm, ConnectionBase::out);
  tester.out() << "Connecting velocity for display ..." << std::endl;
  lux.connect("P-velocity", P.R, P.V, ConnectionBase::out);
  tester.out() << "Connecting E.x and E.y for display ..." << std::endl;
  lux.connect("E-x", EFD.comp(0), ConnectionBase::out);
  lux.connect("E-y", EFD.comp(1), ConnectionBase::out);

  // Wait for everything to be ready to proceed

  tester.out() << "Waiting for ready signal ..." << std::endl;
  lux.ready();
  tester.out() << "Ready complete, moving on." << std::endl;

  // Begin main timestep loop
  for (int it=1; it <= NumSteps; ++it)
  {
    // Advance particle positions
    tester.out() << "Advance particle positions ..." << std::endl;
    P.R = P.R + dt * P.V;

    // Invoke boundary conditions and update particle distribution
    tester.out() << "Synchronize particles ..." << std::endl;
    P.sync(P.R);
   
    // Gather the E field to the particle positions
    tester.out() << "Gather E field ..." << std::endl;
    gather( P.E, EFD, P.R, Particles_t::InterpolatorTag_t() );

    // Advance the particle velocities
    tester.out() << "Advance particle velocities ..." << std::endl;
    P.V = P.V + dt * P.qm * P.E;

    // Update display
    Pooma::blockAndEvaluate();
    tester.out() << "Updating for iters = " << it << std::endl;
    lux.ready();
  }

  // Display the final particle positions, velocities and qm values.
  tester.out() << "PIC2d timestep loop complete!" << std::endl;
  tester.out() << "-----------------------------" << std::endl;

  // Finish up and report results

  int retval = tester.results("Lux Particles PIC2d display test");
  Pooma::finalize();
  return retval;
}

// ACL:rcsinfo
// ----------------------------------------------------------------------
// $RCSfile: lux_test2.cpp,v $   $Author: bfh $
// $Revision: 1.1 $   $Date: 1999/09/17 00:06:59 $
// ----------------------------------------------------------------------
// ACL:rcsinfo