From 7fa30ed23f416cfdeb5446848608079ee838e742 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 Mar 2017 08:04:58 -0800 Subject: Partially fix the paths in the Python wrapper --- Python/classwrapper.i | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) (limited to 'Python/classwrapper.i') diff --git a/Python/classwrapper.i b/Python/classwrapper.i index 7cbf2246a..98f01a88e 100644 --- a/Python/classwrapper.i +++ b/Python/classwrapper.i @@ -24,8 +24,8 @@ %inline %{ std::ifstream & open_ifstream(const char *filename) { std::ifstream *infile = new std::ifstream(filename); - return *infile; - } + return *infile; + } %} %inline %{ @@ -52,7 +52,7 @@ typedef double Real; -%include "../../BoxLib/Src/C_BaseLib/Array.H" +%include "../../amrex/Src/C_BaseLib/Array.H" %extend Array { T& __getitem__ (size_t i) @@ -103,7 +103,7 @@ 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" +//%include "../../amrex/Src/C_BaseLib/IntVect.H" class IntVect { public: @@ -168,20 +168,20 @@ 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/C_BaseLib/Box.H" +//%include "../../amrex/Src/C_BaseLib/FArrayBox.H" +//%include "../../amrex/Src/C_BaseLib/BoxArray.H" +//%include "../../amrex/Src/C_BaseLib/MultiFab.H" //#if (BL_SPACEDIM > 2) %ignore GetDLogA; //#endif -%include "../../BoxLib/Src/C_BaseLib/Geometry.H" +%include "../../amrex/Src/C_BaseLib/Geometry.H" %template(arrayBoxArray) Array; -%include "../../BoxLib/Src/C_ParticleLib/Particles.H" +%include "../../amrex/Src/C_ParticleLib/Particles.H" //%template("WarpXParticleBase") Particle; %template("WarpXParticleContainerBase") ParticleContainer > >; @@ -201,7 +201,7 @@ 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 result(0); self->GetParticleData(result, start_comp, num_comp); @@ -225,6 +225,5 @@ public: } }; -%include "../../BoxLib/Src/C_AmrCoreLib/AmrCore.H" +%include "../../amrex/Src/C_AmrCoreLib/AmrCore.H" %include "../Source/WarpX.H" - -- cgit v1.2.3 From 584a32c8dec1647fa194b5c2db79ef5f61f176c6 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Thu, 9 Mar 2017 14:24:51 -0800 Subject: Converted python interface to AMReX Also added python level time step and incomplete wrap to access the fields. --- Exec/Make.WarpX | 9 +- Python/README | 2 +- Python/classwrapper.i | 157 ++++++--- Python/pywarpx/AMReX.py | 28 ++ Python/pywarpx/BoxLib.py | 28 -- Python/pywarpx/__init__.py | 4 +- Python/pywarpx/timestepper.py | 385 +++++++++++++++++++++++ Python/setup.py | 39 +-- Python/warpxC.i | 5 +- Source/WarpX.H | 36 ++- Source/WarpXEvolve.cpp | 31 ++ tests/Langmuir/langmuir.py | 6 +- tests/plasma_acceleration/plasma_acceleration.py | 6 +- 13 files changed, 616 insertions(+), 120 deletions(-) create mode 100644 Python/pywarpx/AMReX.py delete mode 100644 Python/pywarpx/BoxLib.py create mode 100644 Python/pywarpx/timestepper.py (limited to 'Python/classwrapper.i') diff --git a/Exec/Make.WarpX b/Exec/Make.WarpX index 03a4b8367..7d5c6ebec 100644 --- a/Exec/Make.WarpX +++ b/Exec/Make.WarpX @@ -44,32 +44,31 @@ ifeq ($(USE_PYTHON_MAIN),TRUE) ifeq ($(shell uname),Darwin) SHARED_OPTION = -dynamiclib - MY_LINK_LIBRARIES = $(LDFLAGS) $(libraries) else SHARED_OPTION = -shared - MY_LINK_LIBRARIES = $(LDFLAGS) $(libraries) -lgfortran endif pyinstallswig: libwarpx.a warpxC.i classwrapper.i @echo Making and installing python wrapper ... # touch warpxC.i to force the rebuilding of the swig interface touch warpxC.i - LDFLAGS="$(MY_LINK_LIBRARIES)" CC=$(CXX) CXXFLAGS="$(CXXFLAGS) $(DEFINES)" python setup.py install --force + LDFLAGS="$(LDFLAGS) $(libraries)" CC=$(CXX) CXX=$(CXX) CFLAGS="$(CXXFLAGS) $(CPPFLAGS)" CXXFLAGS="$(CXXFLAGS) $(CPPFLAGS)" DEFINES="$(DEFINES)" python setup.py install --force @echo SUCCESS pybuild: libwarpx.a warpxC.i classwrapper.i @echo Making python wrapper ... - LDFLAGS="$(MY_LINK_LIBRARIES)" CC=$(CXX) CXXFLAGS="$(CXXFLAGS) $(DEFINES)" python setup.py build + LDFLAGS="$(LDFLAGS) $(libraries)" CC=$(CXX) CXX=$(CXX) CFLAGS="$(CXXFLAGS) $(CPPFLAGS)" CXXFLAGS="$(CXXFLAGS) $(CPPFLAGS)" DEFINES="$(DEFINES)" python setup.py build @echo SUCCESS libwarpx.a: $(objForExecs) @echo Making static library $@ ... $(SILENT) $(AR) -crv $@ $^ + $(SILENT) $(RM) AMReX_buildInfo.cpp @echo SUCCESS libwarpx.so: $(objForExecs) @echo Making dynamic library $@ ... - $(SILENT) $(CXX) $(SHARED_OPTION) -fPIC -o $@ $^ $(MY_LINK_LIBRARIES) + $(SILENT) $(CXX) $(SHARED_OPTION) -fPIC -o $@ $^ $(LDFLAGS) $(libraries) $(SILENT) $(RM) AMReX_buildInfo.cpp @echo SUCCESS 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 98f01a88e..8450146ee 100644 --- a/Python/classwrapper.i +++ b/Python/classwrapper.i @@ -1,24 +1,9 @@ %{ -#include -#include -#include -#include -#include - -#include -#include - -#include - -#include -#include -#include -#include -#include - #include +using namespace amrex; + %} %inline %{ @@ -50,11 +35,9 @@ %include %rename(__str__) display; -typedef double Real; +%include "../../amrex/Src/Base/AMReX_Array.H" -%include "../../amrex/Src/C_BaseLib/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; -%template(arrayInt) Array; -%template(arrayGeometry) Array; +// This is needed so that swig knows the amrex::Real is just a float or double +%include "../../amrex/Src/Base/AMReX_REAL.H" + +%template(arrayReal) amrex::Array; +%template(arrayInt) amrex::Array; +%template(arrayGeometry) amrex::Array; -%extend Array { +%extend amrex::Array { 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::size(); + dims[0] = self->std::vector::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 { +%extend amrex::Array { 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::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 "../../amrex/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 "../../amrex/Src/C_BaseLib/Box.H" -//%include "../../amrex/Src/C_BaseLib/FArrayBox.H" -//%include "../../amrex/Src/C_BaseLib/BoxArray.H" -//%include "../../amrex/Src/C_BaseLib/MultiFab.H" +//%include "../../amrex/Src/Base/AMReX_Box.H" +%include "../../amrex/Src/Base/AMReX_BaseFab.H" + +%template(BaseFabReal) amrex::BaseFab; + +%extend amrex::BaseFab { + 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" -//#if (BL_SPACEDIM > 2) +%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 >; +%template(arrayarrayMultifab) amrex::Array > >; + +#if (BL_SPACEDIM > 2) +// GetDLogA is only defined with BL_SPACEDIM <= 2 %ignore GetDLogA; -//#endif +#endif -%include "../../amrex/Src/C_BaseLib/Geometry.H" +%include "../../amrex/Src/Base/AMReX_Geometry.H" -%template(arrayBoxArray) Array; +%template(arrayBoxArray) amrex::Array; -%include "../../amrex/Src/C_ParticleLib/Particles.H" +%include "../../amrex/Src/Particle/AMReX_Particles.H" -//%template("WarpXParticleBase") Particle; -%template("WarpXParticleContainerBase") ParticleContainer > >; +// 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 result(0); + amrex::Array result(0); self->GetParticleLocations(result); npy_intp dims[2] = {BL_SPACEDIM, self->TotalNumberOfParticles()}; PyObject *arr = PyArray_NewFromDescr(&PyArray_Type, @@ -203,7 +255,7 @@ public: return arr; } PyObject * getData(int start_comp, int num_comp) { - Array result(0); + amrex::Array 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 result(0); + amrex::Array result(0); self->GetParticleData(result, 0, num_comp); npy_intp dims[2] = {num_comp, self->TotalNumberOfParticles()}; PyObject *arr = PyArray_NewFromDescr(&PyArray_Type, @@ -225,5 +277,6 @@ public: } }; -%include "../../amrex/Src/C_AmrCoreLib/AmrCore.H" +%include "../../amrex/Src/AmrCore/AMReX_AmrCore.H" %include "../Source/WarpX.H" + diff --git a/Python/pywarpx/AMReX.py b/Python/pywarpx/AMReX.py new file mode 100644 index 000000000..444df3490 --- /dev/null +++ b/Python/pywarpx/AMReX.py @@ -0,0 +1,28 @@ +from .Bucket import Bucket + +from .WarpX import warpx +from .Amr import amr +from .Geometry import geometry +from .Algo import algo +from .Langmuirwave import langmuirwave +from .Interpolation import interpolation +from .Particles import particles + +from . import warpxC + +class AMReX(object): + + def init(self): + argv = [] + argv += warpx.attrlist() + argv += amr.attrlist() + argv += geometry.attrlist() + argv += algo.attrlist() + argv += langmuirwave.attrlist() + argv += interpolation.attrlist() + argv += particles.attrlist() + + warpxC.amrex_init(argv) + + def finalize(self, finalize_mpi=1): + warpxC.amrex_finalize(finalize_mpi) diff --git a/Python/pywarpx/BoxLib.py b/Python/pywarpx/BoxLib.py deleted file mode 100644 index a15119250..000000000 --- a/Python/pywarpx/BoxLib.py +++ /dev/null @@ -1,28 +0,0 @@ -from .Bucket import Bucket - -from .WarpX import warpx -from .Amr import amr -from .Geometry import geometry -from .Algo import algo -from .Langmuirwave import langmuirwave -from .Interpolation import interpolation -from .Particles import particles - -from . import warpxC - -class BoxLib(object): - - def init(self): - argv = [] - argv += warpx.attrlist() - argv += amr.attrlist() - argv += geometry.attrlist() - argv += algo.attrlist() - argv += langmuirwave.attrlist() - argv += interpolation.attrlist() - argv += particles.attrlist() - - warpxC.boxlib_init(argv) - - def finalize(self, finalize_mpi=1): - warpxC.boxlib_finalize(finalize_mpi) diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py index 98c26af36..baec3e39e 100644 --- a/Python/pywarpx/__init__.py +++ b/Python/pywarpx/__init__.py @@ -6,5 +6,7 @@ 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 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 #include #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 %apply (int ARGC, char **ARGV) { (int argc, char *argv[]) } diff --git a/Source/WarpX.H b/Source/WarpX.H index 522e76edc..0b48cbd8b 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -54,6 +54,34 @@ public: static bool use_laser; + amrex::MultiFab * getcurrent (int lev, int direction) {return current[lev][direction].get();} + amrex::MultiFab * getEfield (int lev, int direction) {return Efield[lev][direction].get();} + amrex::MultiFab * getBfield (int lev, int direction) {return Bfield[lev][direction].get();} + + void ComputeDt (); + void MoveWindow (); + + void EvolveE (int lev, amrex::Real dt); + void EvolveB (int lev, amrex::Real dt); + void FillBoundaryE (int lev, bool force); + void FillBoundaryB (int lev, bool force); + void PushParticlesandDepose (int lev, amrex::Real cur_time); + + int getistep (int lev) {return istep[lev];} + void setistep (int lev, int ii) {istep[lev] = ii;} + amrex::Real gett_new (int lev) {return t_new[lev];} + void sett_new (int lev, amrex::Real time) {t_new[lev] = time;} + amrex::Real getdt (int lev) {return dt[lev];} + + int maxStep () {return max_step;} + amrex::Real stopTime () {return stop_time;} + + int checkInt () {return check_int;} + int plotInt () {return plot_int;} + + void WriteCheckPointFile () const; + void WritePlotFile () const; + protected: //! Tagging cells for refinement @@ -88,9 +116,6 @@ private: // Singleton is used when the code is run from python static WarpX* m_instance; - void EvolveE (int lev, amrex::Real dt); - void EvolveB (int lev, amrex::Real dt); - void ReadParameters (); void RegridBaseLevel (); @@ -111,14 +136,9 @@ private: void InitOpenbc (); - void ComputeDt (); - - void MoveWindow (); void InjectPlasma(int num_shift, int dir); void WriteWarpXHeader(const std::string& name) const; - void WriteCheckPointFile() const; - void WritePlotFile () const; void WriteJobInfo (const std::string& dir) const; static void PackPlotDataPtrs (amrex::Array& pmf, const amrex::Array >& data); diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp index e25b2393a..5e5154554 100644 --- a/Source/WarpXEvolve.cpp +++ b/Source/WarpXEvolve.cpp @@ -189,6 +189,37 @@ WarpX::EvolveE (int lev, Real dt) } } +void +WarpX::FillBoundaryE(int lev, bool force) +{ + if (force || WarpX::nox > 1 || WarpX::noy > 1 || WarpX::noz > 1) { + WarpX::FillBoundary(*Efield[lev][0], geom[lev], Ex_nodal_flag); + WarpX::FillBoundary(*Efield[lev][1], geom[lev], Ey_nodal_flag); + WarpX::FillBoundary(*Efield[lev][2], geom[lev], Ez_nodal_flag); + } +} + +void +WarpX::FillBoundaryB(int lev, bool force) +{ + if (force || WarpX::nox > 1 || WarpX::noy > 1 || WarpX::noz > 1) { + WarpX::FillBoundary(*Bfield[lev][0], geom[lev], Bx_nodal_flag); + WarpX::FillBoundary(*Bfield[lev][1], geom[lev], By_nodal_flag); + WarpX::FillBoundary(*Bfield[lev][2], geom[lev], Bz_nodal_flag); + } +} + +void +WarpX::PushParticlesandDepose(int lev, Real cur_time) +{ + // Evolve particles to p^{n+1/2} and x^{n+1} + // Depose current, j^{n+1/2} + mypc->Evolve(lev, + *Efield[lev][0],*Efield[lev][1],*Efield[lev][2], + *Bfield[lev][0],*Bfield[lev][1],*Bfield[lev][2], + *current[lev][0],*current[lev][1],*current[lev][2], cur_time, dt[lev]); +} + void WarpX::ComputeDt () { diff --git a/tests/Langmuir/langmuir.py b/tests/Langmuir/langmuir.py index 8d54a03c1..b0fdd2f8f 100644 --- a/tests/Langmuir/langmuir.py +++ b/tests/Langmuir/langmuir.py @@ -125,8 +125,8 @@ algo.particle_pusher = 0 warpx.cfl = 1.0 # --- Initialize the simulation -boxlib = BoxLib() -boxlib.init() +amrex = AMReX() +amrex.init() warpx.init() set_initial_conditions() @@ -134,4 +134,4 @@ set_initial_conditions() warpx.evolve(max_step) warpx.finalize() -boxlib.finalize() +amrex.finalize() diff --git a/tests/plasma_acceleration/plasma_acceleration.py b/tests/plasma_acceleration/plasma_acceleration.py index 9cbcb1178..f349301e2 100644 --- a/tests/plasma_acceleration/plasma_acceleration.py +++ b/tests/plasma_acceleration/plasma_acceleration.py @@ -223,8 +223,8 @@ warpx.do_plasma_injection = 0 #warpx.injected_plasma_ppc = 4 # --- Initialize the simulation -boxlib = BoxLib() -boxlib.init() +amrex = AMReX() +amrex.init() warpx.init() set_initial_conditions(ncells, domain_min, domain_max) @@ -248,4 +248,4 @@ for i in range(1, max_step + 1): warpx.finalize() -boxlib.finalize() +amrex.finalize() -- cgit v1.2.3