[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re:Re: [pooma-dev] examples/Particles/PIC2d/PIC2d.cpp works
|
From: |
Roman Krylov |
|
Subject: |
Re:Re: [pooma-dev] examples/Particles/PIC2d/PIC2d.cpp works |
|
Date: |
Fri, 09 Jul 2004 18:04:38 +0400 |
|
User-agent: |
Mozilla Thunderbird 0.6 (X11/20040502) |
Perhaps some changes were unnecessary, I was experiencing.
Here they are:
----------
Index: examples/Particles/PIC2d/PIC2d.cpp
===================================================================
RCS file: /home/pooma/Repository/r2/examples/Particles/PIC2d/PIC2d.cpp,v
retrieving revision 1.19
diff -u -r1.19 PIC2d.cpp
--- examples/Particles/PIC2d/PIC2d.cpp 21 Sep 2001 20:27:21 -0000 1.19
+++ examples/Particles/PIC2d/PIC2d.cpp 9 Jul 2004 13:55:04 -0000
@@ -33,16 +33,23 @@
// static electric field using nearest-grid-point interpolation.
//-----------------------------------------------------------------------------
+#include "Field/FieldCentering.h"
+#include "Field/DiffOps/Grad.h"
+#include "Field/DiffOps/Grad.UR.h"
+#include "Particles/InterpolatorNGP.h"
+#include "Particles/Interpolation.h"
#include "Pooma/Particles.h"
#include "Pooma/DynamicArrays.h"
#include "Pooma/Fields.h"
#include "Utilities/Inform.h"
+#include "Pooma/Indices.h"
+
#include <iostream>
#include <stdlib.h>
#include <math.h>
// Traits class for Particles object
-template <class EngineTag, class Centering, class MeshType, class FL,
+template <class EngineTag, class MeshType, class FL,
class InterpolatorTag>
struct PTraits
{
@@ -50,7 +57,7 @@
typedef EngineTag AttributeEngineTag_t;
// The type of particle layout to use
- typedef SpatialLayout<DiscreteGeometry<Centering,MeshType>,FL>
+ typedef SpatialLayout<MeshType,FL>
ParticleLayout_t;
// The type of interpolator to use
@@ -87,6 +94,7 @@
DynamicArray<PointType_t,AttributeEngineTag_t> R;
DynamicArray<PointType_t,AttributeEngineTag_t> V;
DynamicArray<PointType_t,AttributeEngineTag_t> E;
+ DynamicArray<PointType_t,AttributeEngineTag_t> phi;
// for testing purposes of course
DynamicArray<double, AttributeEngineTag_t> qm;
};
@@ -102,24 +110,25 @@
#endif
// Mesh type
-typedef UniformRectilinearMesh<PDim,Cartesian<PDim>,double> Mesh_t;
+typedef UniformRectilinearMesh<PDim,/*,Cartesian<PDim>,*/double> Mesh_t;
// Centering of Field elements on mesh
-typedef Cell Centering_t;
+//typedef CanonicalCentering::CellType Centering_t;
+//static const int Centering_t = CanonicalCentering<PDim>::CellType;
// Geometry type for Fields
-typedef DiscreteGeometry<Centering_t,Mesh_t> Geometry_t;
+//typedef DiscreteGeometry<Centering_t,Mesh_t> Geometry_t;
// Field types
#if POOMA_CHEETAH
-typedef Field< Geometry_t, double,
+typedef Field< Mesh_t, double,
MultiPatch< UniformTag, Remote<Brick> > > DField_t;
-typedef Field< Geometry_t, Vector<PDim,double>,
+typedef Field< Mesh_t, Vector<PDim,double>,
MultiPatch< UniformTag, Remote<Brick> > > VecField_t;
#else
-typedef Field< Geometry_t, double,
+typedef Field< Mesh_t, double,
MultiPatch<UniformTag,Brick> > DField_t;
-typedef Field< Geometry_t, Vector<PDim,double>,
+typedef Field< Mesh_t, Vector<PDim,double>,
MultiPatch<UniformTag,Brick> > VecField_t;
#endif
@@ -131,7 +140,7 @@
typedef NGP InterpolatorTag_t;
// Particle traits class
-typedef PTraits<AttrEngineTag_t,Centering_t,Mesh_t,FLayout_t,
+typedef PTraits<AttrEngineTag_t,Mesh_t,FLayout_t,
InterpolatorTag_t> PTraits_t;
// Type of particle layout
@@ -153,7 +162,7 @@
const double pi = acos(-1.0);
// Maximum value for particle q/m ratio
-const double qmmax = 1.0;
+const double qmmax = 10;//1.0;
// Timestep
const double dt = 1.0;
@@ -169,35 +178,61 @@
out << "-------------------------" << std::endl;
// Create mesh and geometry objects for cell-centered fields.
+ Loc<PDim> blocks(4,4);
+ UniformGridPartition<2> partition(
+ Loc<2>(1, 1),
+ GuardLayers<2>(1), // internal
+ GuardLayers<2>(0)
+ ); // external
Interval<PDim> meshDomain(nx+1,ny+1);
- Mesh_t mesh(meshDomain);
- Geometry_t geometry(mesh);
+ FLayout_t flayout(meshDomain,partition,DistributedTag());
+ Mesh_t mesh(flayout, Vector<2>(0.0), Vector<2>(1.0, 1.0));//(meshDomain);
+ Centering<2> cell =
+ canonicalCentering<2>(CellType, Continuous, AllDim);
+ /*Geometry_t geometry(mesh);*/
// Create a second geometry object that includes a guard layer.
- GuardLayers<PDim> gl(1);
- Geometry_t geometryGL(mesh,gl);
+// GuardLayers<PDim> gl(1);
+// FLayout_t flayoutGL(meshDomain,partition,DistributedTag());
+ /*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,DistributedTag());
- FLayout_t
flayoutGL(geometryGL.physicalDomain(),blocks,gl,DistributedTag());
+// Loc<PDim> blocks(4,4);
+// FLayout_t flayout(mesh.physicalDomain(),blocks,DistributedTag());
+// FLayout_t
flayoutGL(geometryGL.physicalDomain(),blocks,gl,DistributedTag());
// Create and initialize electrostatic potential and electric field.
- DField_t phi(geometryGL,flayoutGL);
- VecField_t EFD(geometry,flayout);
+ DField_t phi(cell,flayout,mesh);
+ VecField_t EFD(cell,flayout,mesh);
// 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);
- Pooma::addAllPeriodicFaceBC(phi, 0.0);
+// Pooma::addAllPeriodicFaceBC(phi, 0.0);
+ phi.addRelation(new Relation0<DField_t,PeriodicFaceBC<PDim>
>(phi,PeriodicFaceBC<PDim>(0)));
+ phi.addRelation(new Relation0<DField_t,PeriodicFaceBC<PDim>
>(phi,PeriodicFaceBC<PDim>(1)));
+ phi.addRelation(new Relation0<DField_t,PeriodicFaceBC<PDim>
>(phi,PeriodicFaceBC<PDim>(2)));
+ phi.addRelation(new Relation0<DField_t,PeriodicFaceBC<PDim>
>(phi,PeriodicFaceBC<PDim>(3)));
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);
+ phi = phi0 * sin(2.0*pi*iota(1,nx).comp(0)/nx)
+ * cos(4.0*pi*iota(1,ny).comp(1)/ny);
+// phi = 100;
+ EFD = -gradVertToCell(phi);
+
+ PrintArray pa;
+ out << "potential: " << std::endl;
+ pa.setDataWidth(10);
+ pa.setScientific(true);
+ pa.print(out.stream(),phi);
+ out << "electric field(to test does grad(phi) work): " << std::endl;
+ pa.setDataWidth(10);
+ pa.setScientific(true);
+ pa.print(out.stream(),EFD);
+
// Create a particle layout object for our use
- PLayout_t layout(geometry,flayout);
+ PLayout_t layout(mesh,flayout);
// Create a Particles object and set periodic boundary conditions
Particles_t P(layout);
@@ -233,7 +268,6 @@
out << "---------------------" << std::endl;
// Display the initial particle positions, velocities and qm values.
- PrintArray pa;
pa.setCarReturn(5);
out << "Initial particle data:" << std::endl;
out << "Particle positions: " << std::endl;
@@ -244,6 +278,11 @@
pa.setDataWidth(10);
pa.setScientific(true);
pa.print(out.stream(),P.V);
+ out << "Field: " << std::endl;
+ pa.setDataWidth(10);
+ pa.setScientific(true);
+ pa.print(out.stream(),P.V);
+
out << "Particle charge-to-mass ratios: " << std::endl;
pa.print(out.stream(),P.qm);
@@ -266,6 +305,7 @@
out << "Advance particle velocities ..." << std::endl;
P.V = P.V + dt * P.qm * P.E;
}
+// gather( P.phi, phi, P.R, Particles_t::InterpolatorTag_t() );//joke :0
// Display the final particle positions, velocities and qm values.
out << "PIC2d timestep loop complete!" << std::endl;
@@ -281,6 +321,11 @@
pa.print(out.stream(),P.V);
out << "Particle charge-to-mass ratios: " << std::endl;
pa.print(out.stream(),P.qm);
+
+ out << "Field: " << std::endl;
+ pa.setDataWidth(10);
+ pa.setScientific(true);
+ pa.print(out.stream(),phi);
// Shut down POOMA and exit
out << "End PIC2d example code." << std::endl;
Index: src/Particles/Interpolation.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/Particles/Interpolation.h,v
retrieving revision 1.7
diff -u -r1.7 Interpolation.h
--- src/Particles/Interpolation.h 26 Oct 2003 12:27:36 -0000 1.7
+++ src/Particles/Interpolation.h 9 Jul 2004 13:55:10 -0000
@@ -50,7 +50,7 @@
//-----------------------------------------------------------------------------
// Includes:
//-----------------------------------------------------------------------------
-
+#include "Evaluator/PatchFunction.h"
//-----------------------------------------------------------------------------
// Forward Declarations:
//-----------------------------------------------------------------------------
Index: src/Particles/InterpolatorNGP.h
===================================================================
RCS file: /home/pooma/Repository/r2/src/Particles/InterpolatorNGP.h,v
retrieving revision 1.13
diff -u -r1.13 InterpolatorNGP.h
--- src/Particles/InterpolatorNGP.h 26 Oct 2003 12:27:36 -0000 1.13
+++ src/Particles/InterpolatorNGP.h 9 Jul 2004 13:55:12 -0000
@@ -160,7 +160,7 @@
// Create a PatchFunction using this functor
PatchFunction< NGPGather<FC,Dim>,
- PatchParticle2<true,false> > patchfun(intfun);
+ PatchParticle2<true,false> > patchfun(intfun);
// Apply the PatchFunction to the attribute using the
// particle position attribute
@@ -444,14 +444,16 @@
Size_t i;
Loc<Dim> indx;
- typedef typename Field_t::Geometry_t Geometry_t;
- const Geometry_t& geom = field_m.geometry();
+ typedef typename Field_t::Mesh_t Mesh_t;
+// typedef typename Field_t::Mesh_t Geometry_t;
+// const Geometry_t& geom = field_m.geometry();
+ const Mesh_t& mesh = field_m.mesh();
for (i=0; i<n; ++i)
{
// Convert the particle position to an index into the Field's
// domain using the Geometry.
- indx = geom.pointIndex(pos(i));
+ indx = mesh.cellContaining(pos(i));
// check we are on the right patch
----------
Richard Guenther wrote:
On Thu, 8 Jul 2004, Roman Krylov wrote:
It's about examples/Particles/PIC2d.
I made it compilable and runnable, though I have no experience in CVS,
and even in RCS so I can't express it that way.
You can get differences compared to the CVS version by issuing
cvs diff -u Interpolation.cpp
f.i., or just
cvs diff -u
for all directories beyond the current one.
I had modified
examples/Particles/PIC2d/PIC2d.cpp,
src/Particles/Interpolation.cpp,
src/Particles/InterpolatorNGP.h,
created 2D spec. of grad(vert2cell only) in src/Field/DiffOps/Grad.h and
Grad.UR.h
In PIC2d.cpp I had to change '<some field>.x().comp(0)' to
'iota(1,nx).comp(0)' for example to make it work and some other subtleties.
Does anybody working on it this time? Perhaps I'm wasting my time if
it's all done already?
I don't know of anyone working with Particles stuff, and I personally are
not very interested in Particles.
But surely it is useful to get Particle - Field interaction back to
working.
Richard.
Roman.
--
Richard Guenther <richard dot guenther at uni-tuebingen dot de>
WWW: http://www.tat.physik.uni-tuebingen.de/~rguenth/