aboutsummaryrefslogtreecommitdiff
path: root/Python
diff options
context:
space:
mode:
Diffstat (limited to 'Python')
-rw-r--r--Python/README2
-rw-r--r--Python/classwrapper.i162
-rw-r--r--Python/pywarpx/AMReX.py (renamed from Python/pywarpx/BoxLib.py)6
-rw-r--r--Python/pywarpx/PGroup.py202
-rw-r--r--Python/pywarpx/__init__.py5
-rw-r--r--Python/pywarpx/timestepper.py385
-rw-r--r--Python/setup.py39
-rw-r--r--Python/warpxC.i5
8 files changed, 727 insertions, 79 deletions
diff --git a/Python/README b/Python/README
index ac70f11fc..13b598d92 100644
--- a/Python/README
+++ b/Python/README
@@ -1,5 +1,5 @@
The initial python wrapped version. This has the minimal capability of allowing specification of the
-input parameters (i.e. the input file) in python, passing that into BoxLib and then running the time steps.
+input parameters (i.e. the input file) in python, passing that into amrex and then running the time steps.
To build it, run "make pyinstallswig" in the Python directory. You will need to have swig installed.
(In the future, the swig generated wrapper could be put into the repo so swig wouldn't be needed.)
diff --git a/Python/classwrapper.i b/Python/classwrapper.i
index 7cbf2246a..8450146ee 100644
--- a/Python/classwrapper.i
+++ b/Python/classwrapper.i
@@ -1,31 +1,16 @@
%{
-#include <iostream>
-#include <ostream>
-#include <fstream>
-#include <sstream>
-#include <stdio.h>
-
-#include <Array.H>
-#include <IntVect.H>
-
-#include <Particles.H>
-
-#include <Box.H>
-#include <FArrayBox.H>
-#include <BoxArray.H>
-#include <MultiFab.H>
-#include <Geometry.H>
-
#include <WarpX.H>
+using namespace amrex;
+
%}
%inline %{
std::ifstream & open_ifstream(const char *filename) {
std::ifstream *infile = new std::ifstream(filename);
- return *infile;
- }
+ return *infile;
+ }
%}
%inline %{
@@ -50,11 +35,9 @@
%include <std_string.i>
%rename(__str__) display;
-typedef double Real;
-
-%include "../../BoxLib/Src/C_BaseLib/Array.H"
+%include "../../amrex/Src/Base/AMReX_Array.H"
-%extend Array {
+%extend amrex::Array {
T& __getitem__ (size_t i)
{
BL_ASSERT(i >= 0);
@@ -73,15 +56,19 @@ typedef double Real;
}
}
-%template(arrayReal) Array<Real>;
-%template(arrayInt) Array<int>;
-%template(arrayGeometry) Array<Geometry>;
+// This is needed so that swig knows the amrex::Real is just a float or double
+%include "../../amrex/Src/Base/AMReX_REAL.H"
-%extend Array<Real> {
+%template(arrayReal) amrex::Array<amrex::Real>;
+%template(arrayInt) amrex::Array<int>;
+%template(arrayGeometry) amrex::Array<amrex::Geometry>;
+
+%extend amrex::Array<amrex::Real> {
PyObject *get_real() {
+ // Get the data as a writeable numpy array, directly accessing the memory.
PyObject *arr = 0;
npy_intp dims[1];
- dims[0] = self->std::vector<Real>::size();
+ dims[0] = self->std::vector<amrex::Real>::size();
arr = PyArray_NewFromDescr(&PyArray_Type,
PyArray_DescrFromType(NPY_DOUBLE), 1, dims, NULL,
self->dataPtr(), NPY_FORTRANORDER|NPY_ARRAY_WRITEABLE, NULL);
@@ -89,8 +76,9 @@ typedef double Real;
return arr;
}
}
-%extend Array<int> {
+%extend amrex::Array<int> {
PyObject *get_int() {
+ // Get the data as a writeable numpy array, directly accessing the memory.
PyObject *arr = 0;
npy_intp dims[1];
dims[0] = self->std::vector<int>::size();
@@ -102,24 +90,31 @@ typedef double Real;
}
}
-// Note that IntVect.H cannot be directly included since swig cannot parse the line setting up "const int* getVect".
-//%include "../../BoxLib/Src/C_BaseLib/IntVect.H"
-class IntVect {
+// This is needed by swig to define the macro D_DECL when including AMReX_IntVect.H
+%include "../../amrex/Src/Base/AMReX_SPACE.H"
+
+// This include can only be done with the modified AMReX_IntVect.H file that hides the ref-qualifiers from swig.
+%include "../../amrex/Src/Base/AMReX_IntVect.H"
+// Save this code, which is needed if AMReX_IntVect.H cannot be included.
+// AMReX_IntVect.H uses ref-qualifiers in the lines setting up getVect that cannot be parsed by swig.
+/*
+class amrex::IntVect {
public:
IntVect(int i, int j, int k);
- IntVect(const IntVect& rhs);
+ IntVect(const amrex::IntVect& rhs);
IntVect& shift(int coord, int s);
// static functions
- static const IntVect& TheZeroVector();
- static const IntVect& TheUnitVector();
- static const IntVect& TheNodeVector();
- static const IntVect& TheCellVector();
+ static const amrex::IntVect& TheZeroVector();
+ static const amrex::IntVect& TheUnitVector();
+ static const amrex::IntVect& TheNodeVector();
+ static const amrex::IntVect& TheCellVector();
};
+*/
-%extend IntVect {
+%extend amrex::IntVect {
//void writeOn(std::ofstream *os){
// *os << *self;
//}
@@ -141,7 +136,7 @@ public:
(*self).setVal(index,val);
}
}
- int __cmp__(const IntVect* other){
+ int __cmp__(const amrex::IntVect* other){
if( (*self) == (*other) ) {
return 0;
}
@@ -168,32 +163,89 @@ public:
//}
}
-//%include "../../BoxLib/Src/C_BaseLib/Box.H"
-//%include "../../BoxLib/Src/C_BaseLib/FArrayBox.H"
-//%include "../../BoxLib/Src/C_BaseLib/BoxArray.H"
-//%include "../../BoxLib/Src/C_BaseLib/MultiFab.H"
+//%include "../../amrex/Src/Base/AMReX_Box.H"
+%include "../../amrex/Src/Base/AMReX_BaseFab.H"
+
+%template(BaseFabReal) amrex::BaseFab<amrex::Real>;
+
+%extend amrex::BaseFab<amrex::Real> {
+ PyObject * get(size_t n=0) {
+ PyObject *arr = 0;
+ npy_intp dims[BL_SPACEDIM];
+ amrex::IntVect size = self->box().size();
+ dims[0] = size[0];
+ dims[1] = size[1];
+ dims[2] = size[2];
+ arr = PyArray_NewFromDescr(&PyArray_Type,
+ PyArray_DescrFromType(NPY_DOUBLE), BL_SPACEDIM, dims, NULL,
+ self->dataPtr(n), NPY_FORTRANORDER|NPY_ARRAY_WRITEABLE, NULL);
+ Py_INCREF(arr);
+ return arr;
+ }
+}
+
+//%include "../../amrex/Src/Base/AMReX_FArrayBox.H"
+//%include "../../amrex/Src/Base/AMReX_BoxArray.H"
+
+%include "../../amrex/Src/Base/AMReX_FabArray.H"
+
+%extend amrex::FabArray {
+ FAB& __getitem__ (size_t i)
+ {
+ BL_ASSERT(i >= 0);
+ BL_ASSERT(i < self->size());
+ return self->get(i);
+ }
+
+ //
+ // Same as above, except acts on const Array's.
+ //
+ const FAB& __getitem__ (size_t i) const
+ {
+ BL_ASSERT(i >= 0);
+ BL_ASSERT(i < self->size());
+ return self->get(i);
+ }
+}
+
+%template(FabArrayFArrayBox) amrex::FabArray< amrex::FArrayBox >;
+
+%include "../../amrex/Src/Base/AMReX_MultiFab.H"
+
+%template(arrayMultifab) amrex::Array< std::unique_ptr<amrex::MultiFab> >;
+%template(arrayarrayMultifab) amrex::Array<amrex::Array< std::unique_ptr<amrex::MultiFab> > >;
-//#if (BL_SPACEDIM > 2)
+#if (BL_SPACEDIM > 2)
+// GetDLogA is only defined with BL_SPACEDIM <= 2
%ignore GetDLogA;
-//#endif
+#endif
-%include "../../BoxLib/Src/C_BaseLib/Geometry.H"
+%include "../../amrex/Src/Base/AMReX_Geometry.H"
-%template(arrayBoxArray) Array<BoxArray>;
+%template(arrayBoxArray) amrex::Array<amrex::BoxArray>;
-%include "../../BoxLib/Src/C_ParticleLib/Particles.H"
+%include "../../amrex/Src/Particle/AMReX_Particles.H"
-//%template("WarpXParticleBase") Particle<PIdx::nattribs,0>;
-%template("WarpXParticleContainerBase") ParticleContainer<PIdx::nattribs,0,std::vector<Particle<PIdx::nattribs,0> > >;
+// Becuase of an apparent problem in swig, the macro around the wrapping of tile_size gives an error during compilation
+%ignore tile_size;
+// Swig doesn't handle the unique_ptr return value causing an error during compilation
%ignore GetChargeDensity;
+// Wrapping WarpXParIter fails since swig doesn't handle the alias SoA properly causing an error during compilation
+%ignore WarpXParIter;
+
+%template("WarpXParticleContainerBase") amrex::ParticleContainer<0,0,PIdx::nattribs>;
+
%include "../Source/ParticleContainer.H"
%include "../Source/WarpXParticleContainer.H"
%extend WarpXParticleContainer {
+ int getNattribs() {
+ return PIdx::nattribs;
+ }
PyObject * getLocations() {
- Array<Real> result(0);
+ amrex::Array<amrex::Real> result(0);
self->GetParticleLocations(result);
npy_intp dims[2] = {BL_SPACEDIM, self->TotalNumberOfParticles()};
PyObject *arr = PyArray_NewFromDescr(&PyArray_Type,
@@ -201,9 +253,9 @@ public:
result.dataPtr(), NPY_ARRAY_F_CONTIGUOUS|NPY_ARRAY_WRITEABLE, NULL);
Py_INCREF(arr);
return arr;
- }
+ }
PyObject * getData(int start_comp, int num_comp) {
- Array<Real> result(0);
+ amrex::Array<amrex::Real> result(0);
self->GetParticleData(result, start_comp, num_comp);
npy_intp dims[2] = {num_comp, self->TotalNumberOfParticles()};
PyObject *arr = PyArray_NewFromDescr(&PyArray_Type,
@@ -214,7 +266,7 @@ public:
}
PyObject * getAllData() {
int num_comp = PIdx::nattribs;
- Array<Real> result(0);
+ amrex::Array<amrex::Real> result(0);
self->GetParticleData(result, 0, num_comp);
npy_intp dims[2] = {num_comp, self->TotalNumberOfParticles()};
PyObject *arr = PyArray_NewFromDescr(&PyArray_Type,
@@ -225,6 +277,6 @@ public:
}
};
-%include "../../BoxLib/Src/C_AmrCoreLib/AmrCore.H"
+%include "../../amrex/Src/AmrCore/AMReX_AmrCore.H"
%include "../Source/WarpX.H"
diff --git a/Python/pywarpx/BoxLib.py b/Python/pywarpx/AMReX.py
index a15119250..444df3490 100644
--- a/Python/pywarpx/BoxLib.py
+++ b/Python/pywarpx/AMReX.py
@@ -10,7 +10,7 @@ from .Particles import particles
from . import warpxC
-class BoxLib(object):
+class AMReX(object):
def init(self):
argv = []
@@ -22,7 +22,7 @@ class BoxLib(object):
argv += interpolation.attrlist()
argv += particles.attrlist()
- warpxC.boxlib_init(argv)
+ warpxC.amrex_init(argv)
def finalize(self, finalize_mpi=1):
- warpxC.boxlib_finalize(finalize_mpi)
+ warpxC.amrex_finalize(finalize_mpi)
diff --git a/Python/pywarpx/PGroup.py b/Python/pywarpx/PGroup.py
new file mode 100644
index 000000000..616141702
--- /dev/null
+++ b/Python/pywarpx/PGroup.py
@@ -0,0 +1,202 @@
+import numpy as np
+from . import WarpX
+from . import warpxC
+
+class PGroup(object):
+ """Implements a class that has the same API as a warp ParticleGroup instance.
+ """
+
+ def __init__(self):
+ self.ns = 1 # Number of species
+ self.npmax = 0 # Size of data arrays
+ self.npid = 0 # number of columns for pid.
+
+ self.gallot()
+
+ def name(self):
+ return 'WarpXParticleGroup'
+
+ def gallot(self):
+ self.lebcancel_pusher = 0 # turns on/off cancellation of E+VxB within V push (logical type)
+ self.lebcancel = 0 # turns on/off cancellation of E+VxB before V push (logical type)
+
+ self.sm = np.zeros(self.ns) # Species mass [kg]
+ self.sq = np.zeros(self.ns) # Species charge [C]
+ self.sw = np.zeros(self.ns) # Species weight, (real particles per simulation particles)
+
+ self.ins = np.ones(self.ns, dtype=int) # Index of first particle in species
+ self.nps = np.zeros(self.ns, dtype=int) # Number of particles in species
+ self.ipmax = np.zeros(self.ns+1, dtype=int) # Max extent within the arrays of each species
+
+ self.sid = np.arange(self.ns, dtype=int) # Global species index for each species
+ self.ndts = np.ones(self.ns, dtype=int) # Stride for time step advance for each species
+ self.ldts = np.ones(self.ns, dtype=int) # (logical type)
+ self.lvdts = np.ones(self.ns, dtype=int) # (logical type)
+ self.iselfb = np.zeros(self.ns, dtype=int) # Group number for particles that are affected by
+ # their own magnetic field, using the 1/gamma**2
+ # approximation. The correction is not applied to
+ # group number -1.
+ self.fselfb = np.zeros(self.ns) # The scaling factor, vz.
+ self.l_maps = np.zeros(self.ns)
+ self.dtscale = np.ones(self.ns) # Scale factor applied to time step size for each
+ # species. Only makes sense in steaday and and
+ # transverse slice modes.
+
+ self.limplicit = np.zeros(self.ns, dtype=int) # Flags implicit particle species (logical type)
+ self.iimplicit = np.full(self.ns, -1, dtype=int) # Group number for implicit particles
+ self.ldoadvance = np.ones(self.ns, dtype=int) # Flags whether particles are time advanced (logical type)
+ self.lboundaries = np.ones(self.ns, dtype=int) # Flags whether boundary conditions need to be applied (logical type)
+ self.lparaxial = np.zeros(self.ns, dtype=int) # Flags to turn on/off paraxial approximation (logical type)
+
+ self.zshift = np.zeros(self.ns)
+ self.gamma_ebcancel_max = np.ones(self.ns) # maximum value allowed for ExB cancellation
+
+ self.gaminv = np.ones(self.npmax) # inverse relativistic gamma factor
+ self._xp = np.zeros(self.npmax) # X-positions of particles [m]
+ self._yp = np.zeros(self.npmax) # Y-positions of particles [m]
+ self._zp = np.zeros(self.npmax) # Z-positions of particles [m]
+ self._uxp = np.zeros(self.npmax) # gamma * X-velocities of particles [m/s]
+ self._uyp = np.zeros(self.npmax) # gamma * Y-velocities of particles [m/s]
+ self._uzp = np.zeros(self.npmax) # gamma * Z-velocities of particles [m/s]
+ self._ex = np.zeros(self.npmax) # Ex of particles [V/m]
+ self._ey = np.zeros(self.npmax) # Ey of particles [V/m]
+ self._ez = np.zeros(self.npmax) # Ez of particles [V/m]
+ self._bx = np.zeros(self.npmax) # Bx of particles [T]
+ self._by = np.zeros(self.npmax) # By of particles [T]
+ self._bz = np.zeros(self.npmax) # Bz of particles [T]
+ self._pid = np.zeros((self.npmax, self.npid)) # Particle ID - used for various purposes
+
+ # --- Temporary fix
+ gchange = gallot
+
+ def allocated(self, name):
+ return (getattr(self, name, None) is not None)
+
+ def addspecies(self):
+ pass
+
+ def _updatelocations(self):
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ xplist = []
+ yplist = []
+ zplist = []
+ for ispecie in range(mypc.nSpecies()):
+ pc = mypc.GetParticleContainer(ispecie)
+ xx = pc.getLocations()
+ xplist.append(xx[0,:])
+ yplist.append(xx[1,:])
+ zplist.append(xx[2,:])
+ self.nps[ispecie] = len(xplist[-1])
+ if ispecie > 0:
+ self.ins[ispecie] = self.ins[ispecie-1] + self.nps[ispecie-1]
+ self.ipmax[ispecie+1] = self.ins[ispecie] + self.nps[ispecie] - 1
+
+ self._xp = np.concatenate(xplist)
+ self._yp = np.concatenate(yplist)
+ self._zp = np.concatenate(zplist)
+ self.npmax = len(self._xp)
+
+ def _updatevelocities(self):
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ uxplist = []
+ uyplist = []
+ uzplist = []
+ for ispecie in range(mypc.nSpecies()):
+ pc = mypc.GetParticleContainer(ispecie)
+ vv = pc.getData(0, 3)
+ uxplist.append(vv[0,:])
+ uyplist.append(vv[1,:])
+ uzplist.append(vv[2,:])
+ self.nps[ispecie] = len(uxplist[-1])
+ if ispecie > 0:
+ self.ins[ispecie] = self.ins[ispecie-1] + self.nps[ispecie-1]
+ self.ipmax[ispecie+1] = self.ins[ispecie] + self.nps[ispecie] - 1
+
+ self._uxp = np.concatenate(uxplist)
+ self._uyp = np.concatenate(uyplist)
+ self._uzp = np.concatenate(uzplist)
+ self.npmax = len(self._xp)
+
+ def _updatepids(self):
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ pidlist = []
+ for ispecie in range(mypc.nSpecies()):
+ pc = mypc.GetParticleContainer(ispecie)
+ self.npid = pc.nAttribs - 3
+ vv = pc.getData(3, self.npid)
+ pidlist.append(vv)
+ self.nps[ispecie] = len(uxplist[-1])
+ if ispecie > 0:
+ self.ins[ispecie] = self.ins[ispecie-1] + self.nps[ispecie-1]
+ self.ipmax[ispecie+1] = self.ins[ispecie] + self.nps[ispecie] - 1
+
+ self._pid = np.concatenate(pidlist.T, axis=0)
+ self.npmax = self._pid.shape[0]
+
+ def getxp(self):
+ self._updatelocations()
+ return self._xp
+ xp = property(getxp)
+
+ def getyp(self):
+ self._updatelocations()
+ return self._yp
+ yp = property(getyp)
+
+ def getzp(self):
+ self._updatelocations()
+ return self._zp
+ zp = property(getzp)
+
+ def getuxp(self):
+ self._updatevelocities()
+ return self._uxp
+ uxp = property(getuxp)
+
+ def getuyp(self):
+ self._updatevelocities()
+ return self._uyp
+ uyp = property(getuyp)
+
+ def getuzp(self):
+ self._updatevelocities()
+ return self._uzp
+ uzp = property(getuzp)
+
+ def getpid(self):
+ self._updatepids()
+ return self._pid
+ pid = property(getpid)
+
+ def getgaminv(self):
+ uxp = self.uxp
+ uyp = self.uyp
+ uzp = self.uzp
+ return sqrt(1. - (uxp**2 + uyp**2 + uzp**2)/warpxC.c**2)
+ gaminv = property(getgaminv)
+
+ def getex(self):
+ return np.zeros(self.npmax)
+ ex = property(getex)
+ def getey(self):
+ return np.zeros(self.npmax)
+ ey = property(getey)
+ def getez(self):
+ return np.zeros(self.npmax)
+ ez = property(getez)
+ def getbx(self):
+ return np.zeros(self.npmax)
+ bx = property(getbx)
+ def getby(self):
+ return np.zeros(self.npmax)
+ by = property(getby)
+ def getbz(self):
+ return np.zeros(self.npmax)
+ bz = property(getbz)
+
diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py
index 98c26af36..4a509419e 100644
--- a/Python/pywarpx/__init__.py
+++ b/Python/pywarpx/__init__.py
@@ -6,5 +6,8 @@ from .Algo import algo
from .Langmuirwave import langmuirwave
from .Interpolation import interpolation
from .Particles import particles
-from .BoxLib import BoxLib
+from .AMReX import AMReX
+
+from .timestepper import TimeStepper
+from .PGroup import PGroup
diff --git a/Python/pywarpx/timestepper.py b/Python/pywarpx/timestepper.py
new file mode 100644
index 000000000..71787c7a4
--- /dev/null
+++ b/Python/pywarpx/timestepper.py
@@ -0,0 +1,385 @@
+import warp
+from warp import top
+from warp import w3d
+
+from . import WarpX
+
+class TimeStepper(object):
+
+
+ def onestep(self):
+ # --- A reference to the instance of the WarpX class
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ self.cur_time = warpx.gett_new(0)
+ self.istep = warpx.getistep(0)
+
+ #if mpi.rank == 0:
+ print "\nSTEP %d starts ..."%(self.istep + 1)
+
+ #if (ParallelDescriptor::NProcs() > 1)
+ # if (okToRegrid(step)) RegridBaseLevel();
+
+ warpx.ComputeDt()
+ dt = warpx.getdt(0)
+
+ # --- Advance level 0 by dt
+ lev = 0
+
+ # --- At the beginning, we have B^{n-1/2} and E^{n}.
+ # --- Particles have p^{n-1/2} and x^{n}.
+ warpx.EvolveB(lev, 0.5*dt) # We now B^{n}
+
+ warpx.FillBoundaryE(lev, False)
+ warpx.FillBoundaryB(lev, False)
+
+ # --- Evolve particles to p^{n+1/2} and x^{n+1}
+ # --- Depose current, j^{n+1/2}
+ warpx.PushParticlesandDepose(lev, self.cur_time)
+
+ mypc.Redistribute() # Redistribute particles
+
+ warpx.EvolveB(lev, 0.5*dt) # We now B^{n+1/2}
+
+ warpx.FillBoundaryB(lev, True)
+
+ warpx.EvolveE(lev, dt) # We now have E^{n+1}
+
+ self.istep += 1
+
+ self.cur_time += dt
+
+ warpx.MoveWindow();
+
+ #if mpi.rank == 0:
+ print "STEP %d ends. TIME = %e DT = %e"%(self.istep, self.cur_time, dt)
+
+ # --- Sync up time
+ for i in range(warpx.finestLevel()+1):
+ warpx.sett_new(i, self.cur_time)
+ warpx.setistep(i, self.istep)
+
+ max_time_reached = ((self.cur_time >= warpx.stopTime() - 1.e-6*dt) or (self.istep >= warpx.maxStep()))
+
+ if warpx.plotInt() > 0 and (self.istep+1)%warpx.plotInt() == 0 or max_time_reached:
+ warpx.WritePlotFile()
+
+ if warpx.checkInt() > 0 and (self.istep+1)%warpx.plotInt() == 0 or max_time_reached:
+ warpx.WriteCheckPointFile()
+
+
+
+
+
+
+
+
+
+
+
+
+
+class TimeStepperFromPICSAR(object):
+
+ def __init__(self, package=None, solver=None, l_debug=False):
+ self.package = package
+ self.solver = solver
+ self.l_debug = l_debug
+
+ def setpackagestepnumber(self, it):
+ if self.package is not None:
+ self.package.setstepnumber(it)
+
+ def step(self, n=1, freq_print=10, lallspecl=0):
+ """
+ This function performs a range of Particle-In-Cell iterations
+
+ Inputs:
+ - n: number of iterations
+ - freq_print: print frequency
+ """
+
+ if (self.l_debug): print("Call step")
+
+ for i in range(n):
+ if(me == 0):
+ if top.it%freq_print == 0:
+ print 'it = %g time = %g'%(top.it, top.time)
+
+ l_first = (lallspecl or (i == 0))
+ l_last = (lallspecl or (i == n-1))
+
+ self.onestep(l_first, l_last)
+
+ if (self.l_debug): print("End step")
+
+ def onestep(self, l_first, l_last):
+ """
+ Perform a single particle-in-cell step
+ """
+
+ if (self.l_debug): print("Call onestep")
+
+ # --- Iteration number
+ self.setpackagestepnumber(top.it)
+
+ # --- call beforestep functions
+ if (self.l_debug): print("Call beforestep functions")
+ warp.callbeforestepfuncs.callfuncsinlist()
+
+ # --- gather fields from grid to particles
+ if (self.l_debug): print("Call Field gathering and particle push")
+
+ # --- push
+ if l_first:
+ if self.package is None:
+ # --- Standard Warp advance
+ for specie in warp.listofallspecies:
+ for pg in specie.flatten(specie.pgroups):
+ for js in range(pg.ns):
+ self.push_velocity_second_part(js, pg)
+ self.push_positions(js, pg)
+ warp.particleboundaries3d(pg, -1, False)
+
+ else:
+ # --- Particle pusher
+ if (self.l_debug): print("Call package particle push")
+ pxr.pxrpush_particles_part2()
+
+ # --- Particle boundary consitions
+ if (self.l_debug): print("Call package particle boundary conditions")
+ pxr.particle_bcs()
+
+ if (self.l_debug): print("Call aliasparticlearrays()")
+ self.aliasparticlearrays()
+
+ else:
+ if self.package is None:
+ # --- Standard Warp advance
+
+ for specie in warp.listofallspecies:
+ for pg in specie.flatten(specie.pgroups):
+ for js in range(pg.ns):
+ self.push_velocity_full(js, pg)
+ self.push_positions(js, pg)
+
+ warp.particleboundaries3d(pg, -1, False)
+
+ else:
+ # --- Particle pusher
+ if (self.l_debug): print("Call package particle pusher")
+ pxr.field_gathering_plus_particle_pusher()
+
+ # --- Particle boundary conditions
+ if (self.l_debug): print("Call package particle boundary conditions")
+ pxr.particle_bcs()
+
+ self.aliasparticlearrays()
+
+ # --- Particle sorting
+ if (self.l_debug): print("Call Particle Sorting")
+ if self.package is not None:
+ # --- This should be a function installed before load rho
+ pxr.particle_sorting_sub()
+
+ # --- call beforeloadrho functions
+ if (self.l_debug): print("Call beforeloadrho functions")
+ warp.beforeloadrho.callfuncsinlist()
+
+ pgroups = []
+ for specie in warp.listofallspecies:
+ pgroups += specie.flatten(specie.pgroups)
+
+ self.pgroups = pgroups
+
+ # --- Call user-defined injection routines
+ if (self.l_debug): print("Call user-defined injection routines")
+ warp.userinjection.callfuncsinlist()
+
+ xgriddiff = w3d.xmmin - pxr.xmin
+ ygriddiff = w3d.ymmin - pxr.ymin
+ zgriddiff = w3d.zmmin - pxr.zmin
+
+ if (xgriddiff != 0 or ygriddiff != 0 or zgriddiff != 0):
+ pxr.pxr_move_sim_boundaries(xgriddiff, ygriddiff, zgriddiff)
+
+
+ if (self.l_debug): print("Call loadrho")
+ self.solver.loadrho(pgroups = pgroups)
+ if (self.l_debug): print("Call loadj")
+ self.solver.loadj(pgroups = pgroups)
+
+ if (self.l_debug): print("Call dosolve")
+ self.solver.dosolve()
+
+ if self.package is None:
+ for specie in warp.listofallspecies:
+ for pg in specie.flatten(specie.pgroups):
+ for js in range(pg.ns):
+ self.fetcheb(js, pg)
+ if l_last:
+ self.push_velocity_first_part(js, pg)
+
+ else:
+ if l_last:
+ if (self.l_debug): print("Call package push particles 1")
+ pxr.pxrpush_particles_part1()
+
+ # --- update time, time counter
+ top.time += top.dt
+ if top.it%top.nhist == 0:
+# zmmnt()
+ minidiag(top.it, top.time, top.lspecial)
+ top.it += 1
+
+ # --- Load balance function should be installed after step
+
+ # --- call afterstep functions
+ if (self.l_debug): print("Call callafterstepfuncs.callfuncsinlist()")
+ warp.callafterstepfuncs.callfuncsinlist()
+
+
+ # --- The following methods use the standard Warp routines
+ def fetcheb(self, js, pg=None):
+ if self.l_verbose:print me, 'enter fetcheb'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ w3d.pgroupfsapi = pg
+ w3d.ipminfsapi = pg.ins[js]
+ w3d.npfsapi = pg.nps[js]
+ pg.ex[il:iu] = 0.
+ pg.ey[il:iu] = 0.
+ pg.ez[il:iu] = 0.
+ pg.bx[il:iu] = 0.
+ pg.by[il:iu] = 0.
+ pg.bz[il:iu] = 0.
+ self.fetche()
+ self.fetchb()
+ w3d.pgroupfsapi = None
+
+ def push_velocity_full(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_velocity_full'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ if pg.lebcancel_pusher:
+ warp.ebcancelpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, 0)
+ else:
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js, pg)
+ # --- push velocity from magnetic field
+ warp.bpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, top.ibpush)
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js, pg)
+
+ if self.l_verbose:print me, 'exit push_velocity_first_part'
+
+ def push_velocity_first_part(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_velocity_first_part'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ if pg.lebcancel_pusher:
+ warp.ebcancelpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, 1)
+ else:
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js,pg)
+ # --- push velocity from magnetic field
+ warp.bpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt, top.ibpush)
+
+ if self.l_verbose:print me, 'exit push_velocity_first_part'
+
+ def push_velocity_second_part(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_velocity_second_part'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ if pg.lebcancel_pusher:
+ warp.ebcancelpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, 2)
+ else:
+ # --- push velocity from magnetic field
+ warp.bpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt, top.ibpush)
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js, pg)
+
+ if self.l_verbose:print me, 'exit push_velocity_second_part'
+
+ def set_gamma(self, js, pg=None):
+ if self.l_verbose:print me, 'enter set_gamma'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ # --- update gamma
+ warp.gammaadv(np, pg.gaminv[il:iu], pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], top.gamadv, top.lrelativ)
+
+ if self.l_verbose:print me, 'exit set_gamma'
+
+ def push_positions(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_positions'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ warp.xpush3d(np, pg.xp[il:iu], pg.yp[il:iu], pg.zp[il:iu],
+ pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.gaminv[il:iu], top.dt)
+
+ if self.l_verbose:print me, 'exit push_positions'
+
+
+
+
+
+
+
+
diff --git a/Python/setup.py b/Python/setup.py
index a9d32134d..3c33a6b3a 100644
--- a/Python/setup.py
+++ b/Python/setup.py
@@ -14,30 +14,33 @@ try:
except AttributeError:
numpy_include = numpy.get_numpy_include()
-boxlib_home = os.environ.get('BOXLIB_HOME', '../../BoxLib')
-boxlib_includes = ['Src/C_BaseLib',
- 'Src/C_ParticleLib',
- 'Src/C_BoundaryLib',
- 'Src/C_AmrCoreLib',
- 'Tools/C_scripts']
-boxlib_includes = [os.path.join(boxlib_home, ii) for ii in boxlib_includes]
-
-include_dirs = [numpy_include, '../Source'] + boxlib_includes
-
-cpp11_flags = ['-std=c++11']
-if platform.system() == "Darwin":
- macosx_deployment_target = platform.mac_ver()[0]
- os.environ['MACOSX_DEPLOYMENT_TARGET'] = macosx_deployment_target
- cpp11_flags.append("-stdlib=libc++")
+amrex_home = os.environ.get('AMREX_HOME', '../../amrex')
+amrex_includes = ['Src/Base',
+ 'Src/Particle',
+ 'Src/Boundary',
+ 'Src/AmrCore',
+ 'Tools/scripts']
+amrex_includes = [os.path.join(amrex_home, ii) for ii in amrex_includes]
+
+include_dirs = [numpy_include, '../Source'] + amrex_includes
+
+definesstring = os.environ.get('DEFINES','')
+defines = definesstring.split(' ')
+
+#cpp11_flags = [] #['-std=c++11']
+#if platform.system() == "Darwin":
+# macosx_deployment_target = platform.mac_ver()[0]
+# os.environ['MACOSX_DEPLOYMENT_TARGET'] = macosx_deployment_target
+# cpp11_flags.append("-stdlib=libc++")
example_module = Extension('pywarpx._warpxC',
- swig_opts=['-c++', '-outdir', 'pywarpx'],
+ swig_opts=['-c++', '-outdir', 'pywarpx'] + defines,
sources=['warpxC.i'],
library_dirs=['.'],
libraries=['warpx'],
include_dirs = include_dirs,
- define_macros = [('BL_USE_MPI','1'), ('BL_SPACEDIM','3'), ('BL_FORT_USE_UNDERSCORE','1'), ('USE_PARTICLES', None)],
- extra_compile_args = cpp11_flags,
+ #define_macros = [('BL_USE_MPI','1'), ('BL_SPACEDIM','3'), ('BL_FORT_USE_UNDERSCORE','1'), ('USE_PARTICLES', None)],
+ #extra_compile_args = cpp11_flags,
)
setup (name = 'pywarpx',
diff --git a/Python/warpxC.i b/Python/warpxC.i
index 0c94829cf..6865a7fb0 100644
--- a/Python/warpxC.i
+++ b/Python/warpxC.i
@@ -1,6 +1,7 @@
%module warpxC
%{
+#include <WarpXConst.H>
#include <WarpXWrappers.h>
#define SWIG_FILE_WITH_INIT
@@ -11,7 +12,9 @@
import_array();
%}
-// For boxlib_init(int argc, char *argv[]);
+%include "../Source/WarpXConst.H"
+
+// For amrex_init(int argc, char *argv[]);
%include <argcargv.i>
%apply (int ARGC, char **ARGV) { (int argc, char *argv[]) }
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-17 21:34:01 +0000