Skip to content

Commit

Permalink
Merge branch 'tuomas-master' int optimize-glue
Browse files Browse the repository at this point in the history
Conflicts:
	fieldsolver/gridGlue.cpp
  • Loading branch information
galfthan committed Apr 24, 2019
2 parents e673413 + a177e3f commit 47fb1ff
Show file tree
Hide file tree
Showing 70 changed files with 1,650 additions and 1,211 deletions.
5 changes: 1 addition & 4 deletions MAKE/Makefile.sisu_gcc
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,4 @@ INC_PROFILE = -I$(LIBRARY_PREFIX)/mpich2/$(MPT_VERSION)/$(CC_BRAND)/$(CC_BRAND_V
INC_EIGEN = -I$(LIBRARY_PREFIX)/eigen/
INC_DCCRG = -I$(LIBRARY_PREFIX)/dccrg_new_neighbours/
INC_VECTORCLASS = -I$(LIBRARY_PREFIX)/vectorclass
#INC_FSGRID = -I$(LIBRARY_PREFIX)/fsgrid
INC_FSGRID = -I/homeappl/home/koskelat/lib/fsgrid/
INC_DCCRG = -I/homeappl/home/koskelat/lib/dccrg/

INC_FSGRID = -I$(LIBRARY_PREFIX)/fsgrid
74 changes: 37 additions & 37 deletions Makefile

Large diffs are not rendered by default.

169 changes: 88 additions & 81 deletions backgroundfield/backgroundfield.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,13 @@
//FieldFunction should be initialized
void setBackgroundField(
FieldFunction& bgFunction,
Real* cellParams,
Real* faceDerivatives,
Real* volumeDerivatives,
FsGrid< std::array<Real, fsgrids::bgbfield::N_BGB>, 2>& BgBGrid,
bool append) {
using namespace CellParams;
using namespace fieldsolver;
using namespace bvolderivatives;

/*if we do not add a new background to the existing one we first put everything to zero*/
if(append==false) {
setBackgroundFieldToZero(BgBGrid);
}

//these are doubles, as the averaging functions copied from Gumics
//use internally doubles. In any case, it should provide more
Expand All @@ -48,19 +48,6 @@ void setBackgroundField(
double dx[3];
unsigned int faceCoord1[3];
unsigned int faceCoord2[3];


start[0] = cellParams[CellParams::XCRD];
start[1] = cellParams[CellParams::YCRD];
start[2] = cellParams[CellParams::ZCRD];

dx[0] = cellParams[CellParams::DX];
dx[1] = cellParams[CellParams::DY];
dx[2] = cellParams[CellParams::DZ];

end[0]=start[0]+dx[0];
end[1]=start[1]+dx[1];
end[2]=start[2]+dx[2];

//the coordinates of the edges face with a normal in the third coordinate direction, stored here to enable looping
faceCoord1[0]=1;
Expand All @@ -69,77 +56,97 @@ void setBackgroundField(
faceCoord2[1]=2;
faceCoord1[2]=0;
faceCoord2[2]=1;

/*if we do not add a new background to the existing one we first put everything to zero*/
if(append==false) {
setBackgroundFieldToZero(cellParams, faceDerivatives, volumeDerivatives);
}

//Face averages
for(unsigned int fComponent=0;fComponent<3;fComponent++){
bgFunction.setDerivative(0);
bgFunction.setComponent((coordinate)fComponent);
cellParams[CellParams::BGBX+fComponent] +=
surfaceAverage(
bgFunction,
(coordinate)fComponent,
accuracy,
start,
dx[faceCoord1[fComponent]],
dx[faceCoord2[fComponent]]
);

//Compute derivatives. Note that we scale by dx[] as the arrays are assumed to contain differences, not true derivatives!
bgFunction.setDerivative(1);
bgFunction.setDerivComponent((coordinate)faceCoord1[fComponent]);
faceDerivatives[fieldsolver::dBGBxdy+2*fComponent] +=
dx[faceCoord1[fComponent]]*
surfaceAverage(bgFunction,(coordinate)fComponent,accuracy,start,dx[faceCoord1[fComponent]],dx[faceCoord2[fComponent]]);
bgFunction.setDerivComponent((coordinate)faceCoord2[fComponent]);
faceDerivatives[fieldsolver::dBGBxdy+1+2*fComponent] +=
dx[faceCoord2[fComponent]]*
surfaceAverage(bgFunction,(coordinate)fComponent,accuracy,start,dx[faceCoord1[fComponent]],dx[faceCoord2[fComponent]]);
}

//Volume averages
for(unsigned int fComponent=0;fComponent<3;fComponent++){
bgFunction.setDerivative(0);
bgFunction.setComponent((coordinate)fComponent);
cellParams[CellParams::BGBXVOL+fComponent] += volumeAverage(bgFunction,accuracy,start,end);

//Compute derivatives. Note that we scale by dx[] as the arrays are assumed to contain differences, not true derivatives!
bgFunction.setDerivative(1);
bgFunction.setDerivComponent((coordinate)faceCoord1[fComponent]);
volumeDerivatives[bvolderivatives::dBGBXVOLdy+2*fComponent] += dx[faceCoord1[fComponent]]*volumeAverage(bgFunction,accuracy,start,end);
bgFunction.setDerivComponent((coordinate)faceCoord2[fComponent]);
volumeDerivatives[bvolderivatives::dBGBXVOLdy+1+2*fComponent] += dx[faceCoord2[fComponent]]*volumeAverage(bgFunction,accuracy,start,end);
auto localSize = BgBGrid.getLocalSize();

#pragma omp parallel for collapse(3)
for (int x = 0; x < localSize[0]; ++x) {
for (int y = 0; y < localSize[1]; ++y) {
for (int z = 0; z < localSize[2]; ++z) {
std::array<double, 3> start3 = BgBGrid.getPhysicalCoords(x, y, z);
start[0] = start3[0];
start[1] = start3[1];
start[2] = start3[2];

dx[0] = BgBGrid.DX;
dx[1] = BgBGrid.DY;
dx[2] = BgBGrid.DZ;

end[0]=start[0]+dx[0];
end[1]=start[1]+dx[1];
end[2]=start[2]+dx[2];

//Face averages
for(uint fComponent=0; fComponent<3; fComponent++){
bgFunction.setDerivative(0);
bgFunction.setComponent((coordinate)fComponent);
BgBGrid.get(x,y,z)->at(fsgrids::bgbfield::BGBX+fComponent) +=
surfaceAverage(bgFunction,
(coordinate)fComponent,
accuracy,
start,
dx[faceCoord1[fComponent]],
dx[faceCoord2[fComponent]]
);

//Compute derivatives. Note that we scale by dx[] as the arrays are assumed to contain differences, not true derivatives!
bgFunction.setDerivative(1);
bgFunction.setDerivComponent((coordinate)faceCoord1[fComponent]);
BgBGrid.get(x,y,z)->at(fsgrids::bgbfield::dBGBxdy+2*fComponent) +=
dx[faceCoord1[fComponent]] *
surfaceAverage(bgFunction,
(coordinate)fComponent,
accuracy,
start,
dx[faceCoord1[fComponent]],
dx[faceCoord2[fComponent]]
);
bgFunction.setDerivComponent((coordinate)faceCoord2[fComponent]);
BgBGrid.get(x,y,z)->at(fsgrids::bgbfield::dBGBxdy+1+2*fComponent) +=
dx[faceCoord2[fComponent]] *
surfaceAverage(bgFunction,
(coordinate)fComponent,
accuracy,
start,
dx[faceCoord1[fComponent]],
dx[faceCoord2[fComponent]]
);
}

//Volume averages
for(unsigned int fComponent=0;fComponent<3;fComponent++){
bgFunction.setDerivative(0);
bgFunction.setComponent((coordinate)fComponent);
BgBGrid.get(x,y,z)->at(fsgrids::bgbfield::BGBXVOL+fComponent) += volumeAverage(bgFunction,accuracy,start,end);

//Compute derivatives. Note that we scale by dx[] as the arrays are assumed to contain differences, not true derivatives!
bgFunction.setDerivative(1);
bgFunction.setDerivComponent((coordinate)faceCoord1[fComponent]);
BgBGrid.get(x,y,z)->at(fsgrids::bgbfield::dBGBXVOLdy+2*fComponent) += dx[faceCoord1[fComponent]] * volumeAverage(bgFunction,accuracy,start,end);
bgFunction.setDerivComponent((coordinate)faceCoord2[fComponent]);
BgBGrid.get(x,y,z)->at(fsgrids::bgbfield::dBGBXVOLdy+1+2*fComponent) += dx[faceCoord2[fComponent]] * volumeAverage(bgFunction,accuracy,start,end);
}
}
}
}

//TODO
//COmpute divergence and curl of volume averaged field and check that both are zero.
}

void setBackgroundFieldToZero(
Real* cellParams,
Real* faceDerivatives,
Real* volumeDerivatives
FsGrid< std::array<Real, fsgrids::bgbfield::N_BGB>, 2>& BgBGrid
) {
using namespace CellParams;
using namespace fieldsolver;
using namespace bvolderivatives;

//Face averages
for(unsigned int fComponent=0;fComponent<3;fComponent++){
cellParams[CellParams::BGBX+fComponent] = 0.0;
faceDerivatives[fieldsolver::dBGBxdy+2*fComponent] = 0.0;
faceDerivatives[fieldsolver::dBGBxdy+1+2*fComponent] = 0.0;
}
auto localSize = BgBGrid.getLocalSize();

//Volume averages
for(unsigned int fComponent=0;fComponent<3;fComponent++){
cellParams[CellParams::BGBXVOL+fComponent] = 0.0;
volumeDerivatives[bvolderivatives::dBGBXVOLdy+2*fComponent] = 0.0;
volumeDerivatives[bvolderivatives::dBGBXVOLdy+1+2*fComponent] =0.0;
#pragma omp parallel for collapse(3)
for (int x = 0; x < localSize[0]; ++x) {
for (int y = 0; y < localSize[1]; ++y) {
for (int z = 0; z < localSize[2]; ++z) {
for (int i = 0; i < fsgrids::bgbfield::N_BGB; ++i) {
BgBGrid.get(x,y,z)->at(i) = 0;
}
}
}
}

}
11 changes: 5 additions & 6 deletions backgroundfield/backgroundfield.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,18 +25,17 @@

#include "fieldfunction.hpp"
#include "../definitions.h"
#include "../common.h"
#include "fsgrid.hpp"

void setBackgroundField(
FieldFunction& bgFunction,
Real* cellParams,
Real* faceDerivatives,
Real* volumeDerivatives,
FsGrid< std::array<Real, fsgrids::bgbfield::N_BGB>, 2>& BgBGrid,
bool append=false
);

void setBackgroundFieldToZero(
Real* cellParams,
Real* faceDerivatives,
Real* volumeDerivatives
FsGrid< std::array<Real, fsgrids::bgbfield::N_BGB>, 2>& BgBGrid
);

#endif
Expand Down
6 changes: 3 additions & 3 deletions fieldsolver/derivatives.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -88,13 +88,13 @@ void calculateDerivatives(
#ifdef DEBUG_SOLVERS
if (leftMoments->at(fsgrids::moments::RHOM) <= 0) {
std::cerr << __FILE__ << ":" << __LINE__
<< (leftMoments->at(fsgrids::moments::RHOM) < 0 ? " Negative" : " Zero") << " density in spatial cell " << leftNbrID
<< (leftMoments->at(fsgrids::moments::RHOM) < 0 ? " Negative" : " Zero") << " density in spatial cell " //<< leftNbrID
<< std::endl;
abort();
}
if (rightMoments->at(fsgrids::moments::RHOM) <= 0) {
if (rghtMoments->at(fsgrids::moments::RHOM) <= 0) {
std::cerr << __FILE__ << ":" << __LINE__
<< (rightMoments->at(fsgrids::moments::RHOM) < 0 ? " Negative" : " Zero") << " density in spatial cell " << rightNbrID
<< (rghtMoments->at(fsgrids::moments::RHOM) < 0 ? " Negative" : " Zero") << " density in spatial cell " //<< rightNbrID
<< std::endl;
abort();
}
Expand Down
Loading

0 comments on commit 47fb1ff

Please sign in to comment.