From 5b54291a04dc68739837d1f37b73faf1bc2dc584 Mon Sep 17 00:00:00 2001 From: Urs Ganse Date: Tue, 19 Jul 2016 11:54:10 +0300 Subject: [PATCH] Re-structure boundary structure in fieldreader for particles. This also affected the fluxfunction (and made that it didn't compile), because it had some hardcoded particleParameters in the field loading code. Now that has been factored out. Also, the fluxfunction tool supports equatorial plane files now (although the support is quite hacky) --- particles/boundaries.h | 256 ++++++++++++++------------- particles/field.h | 36 ++-- particles/particle_post_pusher.cpp | 21 ++- particles/readfields.cpp | 18 +- particles/readfields.h | 13 +- tools/fluxfunction.cpp | 266 ++++++++++++++++++++--------- 6 files changed, 368 insertions(+), 242 deletions(-) diff --git a/particles/boundaries.h b/particles/boundaries.h index 7fde21b31..15c410f24 100644 --- a/particles/boundaries.h +++ b/particles/boundaries.h @@ -2,166 +2,158 @@ #include "particles.h" -class Boundary +struct Boundary { - public: - // Handle a particle's boundary behaviour. - // returns "true" if the particle is still part of the simulation - // afterwards, or "false" if it is to be removed. - virtual bool handleParticle(Particle& p) = 0; - - // Handle cell coordinate in the spatial dimension this boundary - // object cares about (for example: wrap in a periodic direction) - virtual int cellCoordinate(int c) = 0; - - Boundary(int _dimension) : dimension(_dimension) {}; - virtual void setExtent(double _min, double _max, int _cells) { - min=_min; - max=_max; - cells=_cells; - } + // Handle a particle's boundary behaviour. + // returns "true" if the particle is still part of the simulation + // afterwards, or "false" if it is to be removed. + virtual bool handleParticle(Particle& p) = 0; + + // Handle cell coordinate in the spatial dimension this boundary + // object cares about (for example: wrap in a periodic direction) + virtual int cellCoordinate(int c) = 0; + + Boundary(int _dimension) : dimension(_dimension) {}; + virtual void setExtent(double _min, double _max, int _cells) { + min=_min; + max=_max; + cells=_cells; + } - virtual ~Boundary(){}; + virtual ~Boundary(){}; - protected: - // Which spatial dimension to handle - int dimension; + // Which spatial dimension to handle + int dimension; - // Minimum and maximum spatial extents in this dimension - double min, max; + // Minimum and maximum spatial extents in this dimension + double min, max; - // Number of cells in this dimension - int cells; + // Number of cells in this dimension + int cells; }; // Boundary for a spatial dimension that is only 1 cell thick (pseudo-periodic) -class CompactSpatialDimension : public Boundary +struct CompactSpatialDimension : public Boundary { - public: - virtual bool handleParticle(Particle& p) { - // This boundary does not affect particles - return true; - } - virtual int cellCoordinate(int c) { - // Actual cell coordinates in this direction are - // always mapped to 0. - return 0; - } - CompactSpatialDimension(int _dimension) : Boundary(_dimension){}; + virtual bool handleParticle(Particle& p) { + // This boundary does not affect particles + return true; + } + virtual int cellCoordinate(int c) { + // Actual cell coordinates in this direction are + // always mapped to 0. + return 0; + } + CompactSpatialDimension(int _dimension) : Boundary(_dimension){}; }; // Open boundary, which removes particles if they fly out -class OpenBoundary : public Boundary +struct OpenBoundary : public Boundary { - public: - virtual bool handleParticle(Particle& p) { - - // Delete particles that are outside our boundaries. - if(p.x[dimension] <= min || p.x[dimension] >= max) { - return false; - } else { - // Leave all others be. - return true; - } - } + virtual bool handleParticle(Particle& p) { - virtual int cellCoordinate(int c) { - // Cell coordinates are clamped - // TODO: Should this print warnings? - if(c < 0) { - return 0; - } else if(c >= cells) { - return cells-1; - } else { - return c; - } + // Delete particles that are outside our boundaries. + if(p.x[dimension] <= min || p.x[dimension] >= max) { + return false; + } else { + // Leave all others be. + return true; } + } - virtual void setExtent(double _min, double _max, int _cells) { - double dx = (_max-_min)/((double)_cells); - min=_min+2*dx; // 2 Cells border. - max=_max-2*dx; - cells=_cells; + virtual int cellCoordinate(int c) { + // Cell coordinates are clamped + // TODO: Should this print warnings? + if(c < 0) { + return 0; + } else if(c >= cells) { + return cells-1; + } else { + return c; } - OpenBoundary(int _dimension) : Boundary(_dimension){}; + } + + virtual void setExtent(double _min, double _max, int _cells) { + double dx = (_max-_min)/((double)_cells); + min=_min+2*dx; // 2 Cells border. + max=_max-2*dx; + cells=_cells; + } + OpenBoundary(int _dimension) : Boundary(_dimension){}; }; -class ReflectBoundary : public Boundary +struct ReflectBoundary : public Boundary { - public: - virtual bool handleParticle(Particle& p) { - // Particles outside of bounds get their velocities flipped - if(p.x[dimension] <= min || p.x[dimension] >= max) { - p.v *= flip_v; - } - return true; + virtual bool handleParticle(Particle& p) { + // Particles outside of bounds get their velocities flipped + if(p.x[dimension] <= min || p.x[dimension] >= max) { + p.v *= flip_v; } + return true; + } - virtual int cellCoordinate(int c) { - // Cell coordinates are clamped - // TODO: Should this print warnings? - if(c < 0) { - return 0; - } else if(c >= cells) { - return cells-1; - } else { - return c; - } - } - - // Constructor - ReflectBoundary(int _dimension) : Boundary(_dimension) { - double flip[3] = {1.,1.,1.}; - flip[dimension] = -1.; - flip_v.load(flip); - } - virtual void setExtent(double _min, double _max, int _cells) { - double dx = (_max-_min)/((double)_cells); - min=_min+2*dx; // 2 Cells border. - max=_max-2*dx; - cells=_cells; + virtual int cellCoordinate(int c) { + // Cell coordinates are clamped + // TODO: Should this print warnings? + if(c < 0) { + return 0; + } else if(c >= cells) { + return cells-1; + } else { + return c; } - - protected: - // Vector to multiply with in order to flip velocity - // vectors for our dimension - Vec3d flip_v; + } + + // Constructor + ReflectBoundary(int _dimension) : Boundary(_dimension) { + double flip[3] = {1.,1.,1.}; + flip[dimension] = -1.; + flip_v.load(flip); + } + virtual void setExtent(double _min, double _max, int _cells) { + double dx = (_max-_min)/((double)_cells); + min=_min+2*dx; // 2 Cells border. + max=_max-2*dx; + cells=_cells; + } + + // Vector to multiply with in order to flip velocity + // vectors for our dimension + Vec3d flip_v; }; -class PeriodicBoundary : public Boundary +struct PeriodicBoundary : public Boundary { - public: - virtual bool handleParticle(Particle& p) { - if(p.x[dimension] < min) { - p.x += offset_p; - } else if(p.x[dimension] >= max) { - p.x -= offset_p; - } - return true; + virtual bool handleParticle(Particle& p) { + if(p.x[dimension] < min) { + p.x += offset_p; + } else if(p.x[dimension] >= max) { + p.x -= offset_p; } - - virtual int cellCoordinate(int c) { - return c % cells; - } - - // Constructor - PeriodicBoundary(int _dimension) : Boundary(_dimension) { - } - virtual void setExtent(double _min, double _max, int _cells) { - min=_min; - max=_max; - cells=_cells; - - double offset[3] = {0.,0.,0.}; - offset[dimension] = max-min; - offset_p.load(offset); - } - - protected: - // Vector to offset particle positions that leave through - // one boundary with, to come out the other end - Vec3d offset_p; + return true; + } + + virtual int cellCoordinate(int c) { + return c % cells; + } + + // Constructor + PeriodicBoundary(int _dimension) : Boundary(_dimension) { + } + virtual void setExtent(double _min, double _max, int _cells) { + min=_min; + max=_max; + cells=_cells; + + double offset[3] = {0.,0.,0.}; + offset[dimension] = max-min; + offset_p.load(offset); + } + + // Vector to offset particle positions that leave through + // one boundary with, to come out the other end + Vec3d offset_p; }; diff --git a/particles/field.h b/particles/field.h index 06a4793ec..5016ee7ec 100644 --- a/particles/field.h +++ b/particles/field.h @@ -12,36 +12,39 @@ struct Field // Time at which this field is "valid" double time; - // Coordinate boundaries - double min[3]; - double max[3]; - // Mesh spacing double dx[3]; - // Mesh cells - int cells[3]; - - bool periodic[3]; + // Information about spatial dimensions (like extent, boundaries, etc) + Boundary* dimension[3]; // The actual field data std::vector data; + // Constructor (primarily here to make sure boundaries are properly initialized as zero) + Field() { + for(int i=0; i<3; i++) { + dimension[i] = nullptr; + } + } + double* getCellRef(int x, int y, int z) { - if(cells[2] == 1) { - return &(data[4*(y*cells[0]+x)]); + if(dimension[2]->cells == 1) { + // Equatorial plane + return &(data[4*(y*dimension[0]->cells+x)]); } else { - return &(data[4*(z*cells[0]*cells[1] + y*cells[0] + x)]); + // General 3d case + return &(data[4*(z*dimension[0]->cells*dimension[1]->cells + y*dimension[0]->cells + x)]); } } Vec3d getCell(int x, int y, int z) { // Map these cell coordinates using the boundaries - x = ParticleParameters::boundary_behaviour_x->cellCoordinate(x); - y = ParticleParameters::boundary_behaviour_y->cellCoordinate(y); - z = ParticleParameters::boundary_behaviour_z->cellCoordinate(z); + x = dimension[0]->cellCoordinate(x); + y = dimension[1]->cellCoordinate(y); + z = dimension[2]->cellCoordinate(z); double* cell = getCellRef(x,y,z); Vec3d retval; @@ -52,6 +55,7 @@ struct Field // Round-Brace indexing: indexing by physical location, with interpolation Vec3d operator()(Vec3d v) { Vec3d vmin,vdx; + double min[3] = { min[0] = dimension[0]->min, dimension[1]->min, dimension[2]->min}; vmin.load(min); vdx.load(dx); @@ -63,7 +67,7 @@ struct Field truncate_to_int(v).store(index); (v-truncate(v)).store(fract); - if(cells[2] <= 1) { + if(dimension[2]->cells <= 1) { // Equatorial plane Vec3d interp[4]; interp[0] = getCell(index[0],index[1],index[2]); @@ -73,7 +77,7 @@ struct Field return fract[0]*(fract[1]*interp[3]+(1.-fract[1])*interp[1]) + (1.-fract[0])*(fract[1]*interp[2]+(1.-fract[1])*interp[0]); - } else if (cells[1] <= 1) { + } else if (dimension[1]->cells <= 1) { // Polar plane Vec3d interp[4]; diff --git a/particles/particle_post_pusher.cpp b/particles/particle_post_pusher.cpp index a248a3d7e..452d849bf 100644 --- a/particles/particle_post_pusher.cpp +++ b/particles/particle_post_pusher.cpp @@ -38,23 +38,32 @@ int main(int argc, char** argv) { snprintf(filename_buffer,256,filename_pattern.c_str(),input_file_counter-1); readfields(filename_buffer,E[1],B[1],V); E[0]=E[1]; B[0]=B[1]; + E[0].dimension[0] = E[1].dimension[0] = B[0].dimension[0] = B[1].dimension[0] = V.dimension[0] = ParticleParameters::boundary_behaviour_x; + E[0].dimension[1] = E[1].dimension[1] = B[0].dimension[1] = B[1].dimension[1] = V.dimension[1] = ParticleParameters::boundary_behaviour_y; + E[0].dimension[2] = E[1].dimension[2] = B[0].dimension[2] = B[1].dimension[2] = V.dimension[2] = ParticleParameters::boundary_behaviour_z; // Set boundary conditions based on sizes - if(B[0].cells[0] <= 1) { + if(B[0].dimension[0]->cells <= 1) { delete ParticleParameters::boundary_behaviour_x; ParticleParameters::boundary_behaviour_x = createBoundary(0); } - if(B[0].cells[1] <= 1) { + if(B[0].dimension[1]->cells <= 1) { delete ParticleParameters::boundary_behaviour_y; ParticleParameters::boundary_behaviour_y = createBoundary(1); } - if(B[0].cells[2] <= 1) { + if(B[0].dimension[2]->cells <= 1) { delete ParticleParameters::boundary_behaviour_z; ParticleParameters::boundary_behaviour_z = createBoundary(2); } - ParticleParameters::boundary_behaviour_x->setExtent(B[0].min[0], B[0].max[0], B[0].cells[0]); - ParticleParameters::boundary_behaviour_y->setExtent(B[0].min[1], B[0].max[1], B[0].cells[1]); - ParticleParameters::boundary_behaviour_z->setExtent(B[0].min[2], B[0].max[2], B[0].cells[2]); + + // Make sure updated boundary conditions are also correctly known to the fields + E[0].dimension[0] = E[1].dimension[0] = B[0].dimension[0] = B[1].dimension[0] = V.dimension[0] = ParticleParameters::boundary_behaviour_x; + E[0].dimension[1] = E[1].dimension[1] = B[0].dimension[1] = B[1].dimension[1] = V.dimension[1] = ParticleParameters::boundary_behaviour_y; + E[0].dimension[2] = E[1].dimension[2] = B[0].dimension[2] = B[1].dimension[2] = V.dimension[2] = ParticleParameters::boundary_behaviour_z; + + ParticleParameters::boundary_behaviour_x->setExtent(B[0].dimension[0]->min, B[0].dimension[0]->max, B[0].dimension[0]->cells); + ParticleParameters::boundary_behaviour_y->setExtent(B[0].dimension[1]->min, B[0].dimension[1]->max, B[0].dimension[1]->cells); + ParticleParameters::boundary_behaviour_z->setExtent(B[0].dimension[2]->min, B[0].dimension[2]->max, B[0].dimension[2]->cells); /* Init particles */ double dt=ParticleParameters::dt; diff --git a/particles/readfields.cpp b/particles/readfields.cpp index b971cd9ff..54a4fcfee 100644 --- a/particles/readfields.cpp +++ b/particles/readfields.cpp @@ -56,7 +56,7 @@ void debug_output(Field& F, const char* filename) { min[0] = min[1] = min[2] = 99999999999; max[0] = max[1] = max[2] = -99999999999; - for(int i=0; icells*F.dimension[1]->cells; i++) { for(int j=0; j<3; j++) { if(F.data[4*i+j] > max[j]) { max[j] = F.data[4*i+j]; @@ -73,11 +73,11 @@ void debug_output(Field& F, const char* filename) { vmax.load(max); /* Allocate a rgb-pixel array */ - std::vector pixels(4*F.cells[0]*F.cells[1]); + std::vector pixels(4*F.dimension[0]->cells*F.dimension[1]->cells); /* And fill it with colors */ - for(int y=0; ycells; y++) { + for(int x=0; xcells; x++) { /* Rescale the field values to lie between 0..255 */ Vec3d scaled_val = F.getCell(x,y,0); @@ -85,15 +85,15 @@ void debug_output(Field& F, const char* filename) { scaled_val /= (vmax-vmin); scaled_val *= 255.; - pixels[4*(y*F.cells[0] + x)] = (uint8_t) scaled_val[0]; - pixels[4*(y*F.cells[0] + x)+1] = (uint8_t) scaled_val[1]; - pixels[4*(y*F.cells[0] + x)+2] = (uint8_t) scaled_val[2]; - pixels[4*(y*F.cells[0] + x)+3] = 255; // Alpha=1 + pixels[4*(y*F.dimension[0]->cells + x)] = (uint8_t) scaled_val[0]; + pixels[4*(y*F.dimension[0]->cells + x)+1] = (uint8_t) scaled_val[1]; + pixels[4*(y*F.dimension[0]->cells + x)+2] = (uint8_t) scaled_val[2]; + pixels[4*(y*F.dimension[0]->cells + x)+3] = 255; // Alpha=1 } } /* Write it out */ - if(!stbi_write_png(filename, F.cells[0], F.cells[1], 4, pixels.data(), F.cells[0]*4)) { + if(!stbi_write_png(filename, F.dimension[0]->cells, F.dimension[1]->cells, 4, pixels.data(), F.dimension[0]->cells*4)) { std::cerr << "Writing " << filename << " failed: " << strerror(errno) << std::endl; } } diff --git a/particles/readfields.h b/particles/readfields.h index 54ce604d5..b643b7204 100644 --- a/particles/readfields.h +++ b/particles/readfields.h @@ -260,14 +260,21 @@ void readfields(const char* filename, Field& E, Field& B, Field& V, bool doV=tru } } + // Make sure the target fields have boundary data. + if(E.dimension[0] == nullptr || E.dimension[1] == nullptr || E.dimension[2] == nullptr) { + std::cerr << "Warning: Field boundary pointers uninitialized!" << std::endl; + E.dimension[0] = B.dimension[0] = V.dimension[0] = createBoundary(0); + E.dimension[1] = B.dimension[1] = V.dimension[1] = createBoundary(1); + E.dimension[2] = B.dimension[2] = V.dimension[2] = createBoundary(2); + } /* Set field sizes */ for(int i=0; i<3;i++) { /* Volume-centered values -> shift by half a cell in all directions*/ E.dx[i] = B.dx[i] = V.dx[i] = (max[i]-min[i])/cells[i]; double shift = E.dx[i]/2; - E.min[i] = B.min[i] = V.min[i] = min[i]+shift; - E.max[i] = B.max[i] = V.max[i] = max[i]+shift; - E.cells[i] = B.cells[i] = V.cells[i] = cells[i]; + E.dimension[i]->min = B.dimension[i]->min = V.dimension[i]->min = min[i]+shift; + E.dimension[i]->max = B.dimension[i]->max = V.dimension[i]->max = max[i]+shift; + E.dimension[i]->cells = B.dimension[i]->cells = V.dimension[i]->cells = cells[i]; } E.time = B.time = V.time = time; diff --git a/tools/fluxfunction.cpp b/tools/fluxfunction.cpp index 4f4260706..48c448910 100644 --- a/tools/fluxfunction.cpp +++ b/tools/fluxfunction.cpp @@ -12,107 +12,199 @@ using namespace std; -// Calculate fluxfunction by integrating along -z boundary first, -// and then going along z-direction. -std::vector computeFluxUp(Field& B) { - // Create fluxfunction-field to be the same shape as B - std::vector flux(B.cells[0] * B.cells[1] * B.cells[2]); +namespace Polarplane { + // Calculate fluxfunction by integrating along -z boundary first, + // and then going along z-direction. + std::vector computeFluxUp(Field& B) { + // Create fluxfunction-field to be the same shape as B + std::vector flux(B.dimension[0]->cells * B.dimension[1]->cells * B.dimension[2]->cells); - long double tmp_flux=0.; + long double tmp_flux=0.; - // First, fill the z=3 cells - for(int x=B.cells[0]-2; x>0; x--) { - Vec3d bval = B.getCell(x,0,3); + // First, fill the z=3 cells + for(int x=B.dimension[0]->cells-2; x>0; x--) { + Vec3d bval = B.getCell(x,0,3); - tmp_flux -= bval[2] * B.dx[0]; - flux[B.cells[0] * B.cells[1] * 3 + x] = tmp_flux; - } + tmp_flux -= bval[2] * B.dx[0]; + flux[B.dimension[0]->cells * B.dimension[1]->cells * 3 + x] = tmp_flux; + } - // Now, for each row, integrate in z-direction. - for(int x=1; x< B.cells[0]-1; x++) { + // Now, for each row, integrate in z-direction. + for(int x=1; x< B.dimension[0]->cells-1; x++) { - tmp_flux = flux[B.cells[0] * B.cells[1] * 3 + x]; - for(int z=4; z< B.cells[2]; z++) { - Vec3d bval = B.getCell(x,0,z); + tmp_flux = flux[B.dimension[0]->cells * B.dimension[1]->cells * 3 + x]; + for(int z=4; z< B.dimension[2]->cells; z++) { + Vec3d bval = B.getCell(x,0,z); - tmp_flux -= bval[0]*B.dx[2]; - flux[B.cells[0] * B.cells[1] * z + x] = tmp_flux; - } - } + tmp_flux -= bval[0]*B.dx[2]; + flux[B.dimension[0]->cells * B.dimension[1]->cells * z + x] = tmp_flux; + } + } - return flux; -} + return flux; + } -// Calculate fluxfunction by integrating along +z boundary first, -// and then going along negative z-direction. -std::vector computeFluxDown(Field& B) { - // Create fluxfunction-field to be the same shape as B - std::vector flux(B.cells[0] * B.cells[1] * B.cells[2]); + // Calculate fluxfunction by integrating along +z boundary first, + // and then going along negative z-direction. + std::vector computeFluxDown(Field& B) { + // Create fluxfunction-field to be the same shape as B + std::vector flux(B.dimension[0]->cells * B.dimension[1]->cells * B.dimension[2]->cells); - long double tmp_flux=0.; + long double tmp_flux=0.; - // Calculate flux-difference between bottom and top edge - // of +x boundary (so that values are consistent with computeFluxUp) - for(int z=3; zcells-4; z++) { + Vec3d bval = B.getCell(B.dimension[0]->cells-2,0,z); - tmp_flux -= bval[0]*B.dx[2]; - } + tmp_flux -= bval[0]*B.dx[2]; + } - // From there, fill the z=max - 4 cells - for(int x=B.cells[0]-2; x>0; x--) { - Vec3d bval = B.getCell(x,0,B.cells[2]-4); + // First, fill the z=max - 4 cells + for(int x=B.dimension[0]->cells-2; x>0; x--) { + Vec3d bval = B.getCell(x,0,B.dimension[2]->cells-4); - tmp_flux -= bval[2] * B.dx[0]; - flux[B.cells[0] * B.cells[1] * (B.cells[2] - 4) + x] = tmp_flux; - } + tmp_flux -= bval[2] * B.dx[0]; + flux[B.dimension[0]->cells * B.dimension[1]->cells * (B.dimension[2]->cells - 4) + x] = tmp_flux; + } - // Now, for each row, integrate in -z-direction. - for(int x=1; x< B.cells[0]-1; x++) { + // Now, for each row, integrate in -z-direction. + for(int x=1; x< B.dimension[0]->cells-1; x++) { - tmp_flux = flux[B.cells[0] * B.cells[1] * (B.cells[2] - 4) + x]; - for(int z=B.cells[2]-5; z > 0; z--) { - Vec3d bval = B.getCell(x,0,z); + tmp_flux = flux[B.dimension[0]->cells * B.dimension[1]->cells * (B.dimension[2]->cells - 4) + x]; + for(int z=B.dimension[2]->cells-5; z > 0; z--) { + Vec3d bval = B.getCell(x,0,z); - tmp_flux += bval[0] * B.dx[2]; - flux[B.cells[0] * B.cells[1] * z + x] = tmp_flux; - } - } + tmp_flux += bval[0] * B.dx[2]; + flux[B.dimension[0]->cells * B.dimension[1]->cells * z + x] = tmp_flux; + } + } - return flux; -} + return flux; + } -// Calculate fluxfunction by integrating along -x from the right boundary -std::vector computeFluxLeft(Field& B) { - // Create fluxfunction-field to be the same shape as B - std::vector flux(B.cells[0] * B.cells[1] * B.cells[2]); + // Calculate fluxfunction by integrating along -x from the right boundary + std::vector computeFluxLeft(Field& B) { + // Create fluxfunction-field to be the same shape as B + std::vector flux(B.dimension[0]->cells * B.dimension[1]->cells * B.dimension[2]->cells); - long double tmp_flux=0.; - long double bottom_right_flux=0.; + long double tmp_flux=0.; + long double bottom_right_flux=0.; - // First calculate flux difference to bottom right corner - // Now, for each row, integrate in -z-direction. - for(int z=0; z < B.cells[2]; z++) { - Vec3d bval = B.getCell(B.cells[0]-1,0,z); - bottom_right_flux -= bval[0] * B.dx[2]; + // First calculate flux difference to bottom right corner + // Now, for each row, integrate in -z-direction. + for(int z=0; z < B.dimension[2]->cells; z++) { + Vec3d bval = B.getCell(B.dimension[0]->cells-1,0,z); + bottom_right_flux -= bval[0] * B.dx[2]; - tmp_flux = bottom_right_flux; - for(int x=B.cells[0]-1; x>0; x--) { + tmp_flux = bottom_right_flux; + for(int x=B.dimension[0]->cells-1; x>0; x--) { - bval = B.getCell(x,0,z); + bval = B.getCell(x,0,z); - tmp_flux -= bval[2] * B.dx[0]; - flux[B.cells[0] * B.cells[1] * z + x] = tmp_flux; - } - } + tmp_flux -= bval[2] * B.dx[0]; + flux[B.dimension[0]->cells * B.dimension[1]->cells * z + x] = tmp_flux; + } + } - return flux; + return flux; + } } +namespace Equatorialplane { + // Calculate fluxfunction by integrating along -y boundary first, + // and then going along y-direction. + std::vector computeFluxUp(Field& B) { + // Create fluxfunction-field to be the same shape as B + std::vector flux(B.dimension[0]->cells * B.dimension[1]->cells * B.dimension[2]->cells); + + long double tmp_flux=0.; + + // First, fill the y=3 cells + for(int x=B.dimension[0]->cells-2; x>0; x--) { + Vec3d bval = B.getCell(x,3,0); + + tmp_flux -= bval[1] * B.dx[0]; + flux[B.dimension[0]->cells * 3 + x] = tmp_flux; + } + + // Now, for each row, integrate in y-direction. + for(int x=1; x< B.dimension[0]->cells-1; x++) { + + tmp_flux = flux[B.dimension[0]->cells * 3 + x]; + for(int y=4; y< B.dimension[1]->cells; y++) { + Vec3d bval = B.getCell(x,y,0); + + tmp_flux -= bval[0]*B.dx[1]; + flux[B.dimension[0]->cells * y + x] = tmp_flux; + } + } + + return flux; + } + + + + // Calculate fluxfunction by integrating along +y boundary first, + // and then going along negative y-direction. + std::vector computeFluxDown(Field& B) { + // Create fluxfunction-field to be the same shape as B + std::vector flux(B.dimension[0]->cells * B.dimension[1]->cells * B.dimension[2]->cells); + + long double tmp_flux=0.; + + // First, fill the y=max - 4 cells + for(int x=B.dimension[0]->cells-2; x>0; x--) { + Vec3d bval = B.getCell(x,B.dimension[1]->cells-4,0); + + tmp_flux -= bval[1] * B.dx[0]; + flux[B.dimension[0]->cells * (B.dimension[1]->cells - 4) + x] = tmp_flux; + } + + // Now, for each row, integrate in -y-direction. + for(int x=1; x< B.dimension[0]->cells-1; x++) { + + tmp_flux = flux[B.dimension[0]->cells * (B.dimension[1]->cells - 4) + x]; + for(int y=B.dimension[1]->cells-5; y > 0; y--) { + Vec3d bval = B.getCell(x,y,0); + + tmp_flux += bval[0] * B.dx[1]; + flux[B.dimension[0]->cells * y + x] = tmp_flux; + } + } + + return flux; + } + + + + // Calculate fluxfunction by integrating along -x from the right boundary + std::vector computeFluxLeft(Field& B) { + // Create fluxfunction-field to be the same shape as B + std::vector flux(B.dimension[0]->cells * B.dimension[1]->cells * B.dimension[2]->cells); + + long double tmp_flux=0.; + + // Now, for each row, integrate in -y-direction. + for(int y=0; y < B.dimension[1]->cells; y++) { + tmp_flux = 0; + for(int x=B.dimension[0]->cells-1; x>0; x--) { + + Vec3d bval = B.getCell(x,y,0); + + tmp_flux -= bval[1] * B.dx[0]; + flux[B.dimension[0]->cells * y + x] = tmp_flux; + } + } + + return flux; + } + +} // Get a median of 3 values (branch-free!) @@ -137,11 +229,33 @@ int main(int argc, char** argv) { // TODO: Don't uselessly read E, we really only care about B. Field E,B,V; + bool eqPlane = false; readfields(inFile.c_str(),E,B,V,false); + // Make sure we are working with a 2D simulation here. + if(B.dimension[0]->cells > 1 && B.dimension[1]->cells > 1 && B.dimension[2]->cells > 1) { + cerr << "This is a 3D simulation output. Flux function calculation only makes sense for 2D data." + << endl; + exit(1); + } + // Check if we are in x-y plane + if(B.dimension[1]->cells > 1) { + eqPlane = true; + } + cerr << "File read, calculating flux function..." << endl; - std::vector fluxUp(computeFluxUp(B)), fluxDown(computeFluxDown(B)), fluxLeft(computeFluxLeft(B)); + + std::vector fluxUp, fluxDown, fluxLeft; + if(eqPlane) { + fluxUp = Equatorialplane::computeFluxUp(B); + fluxDown = Equatorialplane::computeFluxDown(B); + fluxLeft = Equatorialplane::computeFluxLeft(B); + } else { + fluxUp = Polarplane::computeFluxUp(B); + fluxDown = Polarplane::computeFluxDown(B); + fluxLeft = Polarplane::computeFluxLeft(B); + } for(unsigned int i=0; icells*B.dimension[1]->cells*B.dimension[2]->cells*sizeof(double); // Write binary blob for(ssize_t remain=size; remain > 0; ) { @@ -177,10 +291,10 @@ int main(int argc, char** argv) { } fprintf(f, "TIME: %lf\n", B.time); fprintf(f, "DATA_FILE: %s\n", outFile.c_str()); - fprintf(f, "DATA_SIZE: %i %i %i\n", B.cells[0], B.cells[1], B.cells[2]); + fprintf(f, "DATA_SIZE: %i %i %i\n", B.dimension[0]->cells, B.dimension[1]->cells, B.dimension[2]->cells); fprintf(f, "DATA_FORMAT: DOUBLE\nVARIABLE: fluxfunction\nDATA_ENDIAN: LITTLE\nCENTERING: zonal\n"); - fprintf(f, "BRICK_ORIGIN: %lf %lf %lf\n", B.min[0], B.min[1], B.min[2]); - fprintf(f, "BRICK_SIZE: %lf %lf %lf\n", B.max[0] - B.min[0], B.max[1] - B.min[1], B.max[2] - B.min[2]); + fprintf(f, "BRICK_ORIGIN: %lf %lf %lf\n", B.dimension[0]->min, B.dimension[1]->min, B.dimension[2]->min); + fprintf(f, "BRICK_SIZE: %lf %lf %lf\n", B.dimension[0]->max - B.dimension[0]->min, B.dimension[1]->max - B.dimension[1]->min, B.dimension[2]->max - B.dimension[2]->min); fprintf(f, "DATA_COMPONENTS: 1\n"); fclose(f);