From baa1b5cfd8cf383f811fb9b17e04696812060d30 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 1 Mar 2010 14:56:27 +0100 Subject: [PATCH] Operator as 3D array, Engine and Excitation --- FDTD/cartoperator.cpp | 89 +++++++++++++++++++++++------------------- FDTD/engine.cpp | 91 ++++++++++++++++++++++++++++++++++++++++++- FDTD/engine.h | 13 ++++++- FDTD/operator.cpp | 58 ++++++++++++++++++++++----- FDTD/operator.h | 21 +++++++--- main.cpp | 22 +++++++++++ openEMS.pro | 6 ++- tools/array_ops.cpp | 58 +++++++++++++++++++++++++++ tools/array_ops.h | 12 ++++++ tools/constants.h | 1 + 10 files changed, 310 insertions(+), 61 deletions(-) create mode 100644 tools/array_ops.cpp create mode 100644 tools/array_ops.h diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp index b5fe0f0..ec919ec 100644 --- a/FDTD/cartoperator.cpp +++ b/FDTD/cartoperator.cpp @@ -1,15 +1,14 @@ #include "cartoperator.h" +#include "tools/array_ops.h" CartOperator::CartOperator() { Init(); - Operator::Init(); } CartOperator::~CartOperator() { Reset(); - Operator::Reset(); } void CartOperator::Init() @@ -24,7 +23,6 @@ void CartOperator::Init() EC_L[n]=NULL; EC_R[n]=NULL; } - Operator::Init(); } @@ -39,6 +37,7 @@ void CartOperator::Reset() delete[] EC_L[n]; delete[] EC_R[n]; } + Operator::Reset(); Init(); } @@ -71,25 +70,34 @@ int CartOperator::CalcECOperator() CalcTimestep(); + Delete_N_3DArray(vv,numLines); + Delete_N_3DArray(vi,numLines); + Delete_N_3DArray(iv,numLines); + Delete_N_3DArray(ii,numLines); + vv = Create_N_3DArray(numLines); + vi = Create_N_3DArray(numLines); + iv = Create_N_3DArray(numLines); + ii = Create_N_3DArray(numLines); + + unsigned int i=0; + unsigned int pos[3]; + for (int n=0;n<3;++n) { - delete[] vv[n]; - vv[n] = new FDTD_FLOAT[MainOp->GetSize()]; - delete[] vi[n]; - vi[n] = new FDTD_FLOAT[MainOp->GetSize()]; - delete[] iv[n]; - iv[n] = new FDTD_FLOAT[MainOp->GetSize()]; - delete[] ii[n]; - ii[n] = new FDTD_FLOAT[MainOp->GetSize()]; - - for (unsigned int i=0;iGetSize();++i) + for (pos[0]=0;pos[0]SetPos(pos[0],pos[1],pos[2]); + vv[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); + vi[n][pos[0]][pos[1]][pos[2]] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); - ii[n][i] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); - iv[n][i] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); -// cerr << iv[n][i] << endl; + ii[n][pos[0]][pos[1]][pos[2]] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + iv[n][pos[0]][pos[1]][pos[2]] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + } + } } } @@ -125,14 +133,11 @@ void CartOperator::ApplyElectricBC(bool* dirs) for (pos[nPP]=0;pos[nPP]SetPos(pos[0],pos[1],pos[2]); - vv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n]; - vi[n][ipos] *= (FDTD_FLOAT)!dirs[2*n]; - + vv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; + vi[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; pos[n]=numLines[n]-1; - ipos=MainOp->SetPos(pos[0],pos[1],pos[2]); - vv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1]; - vi[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1]; + vv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; + vi[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; } } } @@ -152,14 +157,12 @@ void CartOperator::ApplyMagneticBC(bool* dirs) for (pos[nPP]=0;pos[nPP]SetPos(pos[0],pos[1],pos[2]); - ii[n][ipos] *= (FDTD_FLOAT)!dirs[2*n]; - iv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n]; + ii[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; + iv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n]; pos[n]=numLines[n]-2; - ipos=MainOp->SetPos(pos[0],pos[1],pos[2]); - ii[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1]; - iv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1]; + ii[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; + iv[n][pos[0]][pos[1]][pos[2]] *= (FDTD_FLOAT)!dirs[2*n+1]; } } } @@ -407,9 +410,10 @@ double CartOperator::CalcTimestep() bool CartOperator::CalcEFieldExcitation() { - vector vIndex; + if (dT==0) return false; + vector vIndex[3]; vector vExcit[3]; - vector vDelay; + vector vDelay; unsigned int ipos; unsigned int pos[3]; double coord[3]; @@ -420,7 +424,6 @@ bool CartOperator::CalcEFieldExcitation() { for (pos[0]=0;pos[0]SetPos(pos[0],pos[1],pos[2]); coord[0] = discLines[0][pos[0]]; coord[1] = discLines[1][pos[1]]; coord[2] = discLines[2][pos[2]]; @@ -430,9 +433,10 @@ bool CartOperator::CalcEFieldExcitation() CSPropElectrode* elec = prop->ToElectrode(); if (elec->GetExcitType()==0) { - vIndex.push_back(ipos); + vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); for (int n=0;n<3;++n) { + vIndex[n].push_back(pos[n]); double delta=MainOp->GetIndexDelta(n,pos[n])*gridDelta; vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta); } @@ -441,13 +445,18 @@ bool CartOperator::CalcEFieldExcitation() } } } - E_Ex_Count = vIndex.size(); - delete[] E_Ex_index; - E_Ex_index = new unsigned int[E_Ex_Count]; + E_Ex_Count = vIndex[0].size(); + for (int n=0;n<3;++n) + { + delete[] E_Ex_index[n]; + E_Ex_index[n] = new unsigned int[E_Ex_Count]; + for (unsigned int i=0;inumLines); + curr = Create_N_3DArray(Op->numLines); +} + +void Engine::Reset() +{ + Delete_N_3DArray(volt,Op->numLines); + volt=NULL; + Delete_N_3DArray(curr,Op->numLines); + curr=NULL; +} + +bool Engine::IterateTS(unsigned int iterTS) +{ + unsigned int pos[3]; + int exc_pos; + + for (unsigned int iter=0;iternumLines[2]-1;++pos[2]) + { + for (pos[1]=1;pos[1]numLines[1]-1;++pos[1]) + { + for (pos[0]=1;pos[0]numLines[0]-1;++pos[0]) + { + //do the updates here + //for x + volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]]; + volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-1][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-1]); + + //for y + volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]]; + volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-1] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-1][pos[1]][pos[2]]); + + //for x + volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]]; + volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-1][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-1][pos[2]]); + } + } + } + + //soft voltage excitation here (E-field excite) + for (unsigned int n=0;nE_Ex_Count;++n) + { + exc_pos = (int)numTS - (int)Op->E_Ex_delay[n]; + exc_pos*= (exc_pos>0 && exc_posExciteLength); + volt[0][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[0][n]*Op->ExciteSignal[exc_pos]; + volt[1][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[1][n]*Op->ExciteSignal[exc_pos]; + volt[2][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[2][n]*Op->ExciteSignal[exc_pos]; + } + + //current updates + for (pos[0]=0;pos[0]numLines[0]-1;++pos[0]) + { + for (pos[1]=0;pos[1]numLines[1]-1;++pos[1]) + { + for (pos[2]=0;pos[2]numLines[2]-1;++pos[2]) + { + //do the updates here + //for x + curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]]; + curr[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]+1][pos[2]] - volt[2][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]+1] + volt[1][pos[0]][pos[1]][pos[2]]); + + //for y + curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]]; + curr[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]+1] - volt[0][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]+1][pos[1]][pos[2]] + volt[2][pos[0]][pos[1]][pos[2]]); + + //for x + curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]]; + curr[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]+1][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]+1][pos[2]] + volt[0][pos[0]][pos[1]][pos[2]]); + } + } + } + + //soft current excitation here (H-field excite) + ++numTS; + } + return true; } diff --git a/FDTD/engine.h b/FDTD/engine.h index 07bca4a..57db067 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -1,17 +1,26 @@ #ifndef ENGINE_H #define ENGINE_H +#include "operator.h" + class Engine { public: - Engine(); + Engine(Operator* op); virtual ~Engine(); + virtual void Init(); + virtual void Reset(); + //!Iterate a number of timesteps - bool IterateTS(unsigned int numTS); + virtual bool IterateTS(unsigned int iterTS); protected: + Operator* Op; + FDTD_FLOAT**** volt; + FDTD_FLOAT**** curr; + unsigned int numTS; }; #endif // ENGINE_H diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index a138aad..5cd6b3e 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -1,4 +1,5 @@ #include "operator.h" +#include "tools/array_ops.h" Operator::Operator() { @@ -14,30 +15,33 @@ void Operator::Init() { CSX = NULL; - E_Ex_index = NULL; + ExciteSignal = NULL; E_Ex_delay = NULL; + vv=NULL; + vi=NULL; + iv=NULL; + ii=NULL; for (int n=0;n<3;++n) { discLines[n]=NULL; - vv[n]=NULL; - vi[n]=NULL; - iv[n]=NULL; - ii[n]=NULL; E_Ex_amp[n]=NULL; + E_Ex_index[n]=NULL; } } void Operator::Reset() { - delete[] E_Ex_index; + delete[] ExciteSignal; delete[] E_Ex_delay; + Delete_N_3DArray(vv,numLines); + Delete_N_3DArray(vi,numLines); + Delete_N_3DArray(iv,numLines); + Delete_N_3DArray(ii,numLines); for (int n=0;n<3;++n) { - delete[] vv[n]; - delete[] vi[n]; - delete[] iv[n]; - delete[] ii[n]; + delete[] discLines[n]; delete[] E_Ex_amp[n]; + delete[] E_Ex_index[n]; } Operator::Init(); } @@ -49,3 +53,37 @@ void Operator::SetGeometryCSX(ContinuousStructure* geo) Reset(); CSX = geo; } + +double Operator::GetNumberCells() +{ + if (numLines) + return (numLines[0]-1)*(numLines[1]-1)*(numLines[2]-1); + return 0; +} + +void Operator::ShowSize() +{ + unsigned int OpSize = 12*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT); + unsigned int FieldSize = 6*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT); + double MBdiff = 1024*1024; + + cout << "FDTD Operator Size:" << endl; + cout << "Size of Operator in Byte : " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MB) " << endl; + cout << "Size of Field-Data in Byte: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MB) " << endl; +} + +void Operator::CalcGaussianPulsExcitation(double f0, double fc) +{ + if (dT==0) return; + + ExciteLength = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT); + cerr << "Operator::CalcGaussianPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl; + delete[] ExciteSignal; + ExciteSignal = new FDTD_FLOAT[ExciteLength+1]; + ExciteSignal[0]=0.0; + for (unsigned int n=1;n #include #include "FDTD/cartoperator.h" +#include "FDTD/engine.h" + #include "ContinuousStructure.h" using namespace std; @@ -21,10 +23,29 @@ int main(int argc, char *argv[]) cop.CalcECOperator(); + cop.CalcGaussianPulsExcitation(1e9,1e9); + time_t OpDoneTime=time(NULL); + cop.ShowSize(); + cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl; + Engine eng(&cop); + + time_t currTime = time(NULL); + + unsigned int NrIter = 500; + + eng.IterateTS(NrIter); + + time_t prevTime = currTime; + currTime = time(NULL); + + double t_diff = difftime(currTime,prevTime); + + cerr << "Time for " << NrIter << " iterations with " << cop.GetNumberCells() << " cells : " << t_diff << " sec" << endl; + cerr << "Speed (MCells/s): " << (double)cop.GetNumberCells()*(double)NrIter/t_diff/1e6 << endl; return 0; } @@ -53,6 +74,7 @@ void BuildMSL(ContinuousStructure &CSX) CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); elec->SetExcitation(1,1); elec->SetExcitType(0); +// elec->SetDelay(2.0e-9); CSX.AddProperty(elec); box = new CSPrimBox(CSX.GetParameterSet(),elec); diff --git a/openEMS.pro b/openEMS.pro index bbd644a..c50de51 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -19,10 +19,12 @@ SOURCES += main.cpp \ tools/ErrorMsg.cpp \ tools/AdrOp.cpp \ FDTD/engine.cpp \ - FDTD/operator.cpp + FDTD/operator.cpp \ + tools/array_ops.cpp HEADERS += FDTD/cartoperator.h \ tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ FDTD/engine.h \ - FDTD/operator.h + FDTD/operator.h \ + tools/array_ops.h diff --git a/tools/array_ops.cpp b/tools/array_ops.cpp new file mode 100644 index 0000000..637708a --- /dev/null +++ b/tools/array_ops.cpp @@ -0,0 +1,58 @@ +#include "array_ops.h" + +FDTD_FLOAT*** Create3DArray(unsigned int* numLines) +{ + FDTD_FLOAT*** array; + unsigned int pos[3]; + array = new FDTD_FLOAT**[numLines[0]]; + for (pos[0]=0;pos[0]