Operator as 3D array, Engine and Excitation

pull/1/head
Thorsten Liebig 2010-03-01 14:56:27 +01:00
parent 50e8ddaf0f
commit baa1b5cfd8
10 changed files with 310 additions and 61 deletions

View File

@ -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;i<MainOp->GetSize();++i)
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
{
vv[n][i] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
vi[n][i] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
{
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
{
i = MainOp->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]<numLines[nPP];++pos[nPP])
{
pos[n]=0;
ipos=MainOp->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]<numLines[nPP];++pos[nPP])
{
pos[n]=0;
ipos=MainOp->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<unsigned int> vIndex;
if (dT==0) return false;
vector<unsigned int> vIndex[3];
vector<FDTD_FLOAT> vExcit[3];
vector<FDTD_FLOAT> vDelay;
vector<unsigned int> vDelay;
unsigned int ipos;
unsigned int pos[3];
double coord[3];
@ -420,7 +424,6 @@ bool CartOperator::CalcEFieldExcitation()
{
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
{
ipos = MainOp->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];
delete[] E_Ex_delay;
E_Ex_delay = new FDTD_FLOAT[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;i<E_Ex_Count;++i)
E_Ex_delay[i]=vIndex.at(i);
E_Ex_index[n][i]=vIndex[n].at(i);
}
delete[] E_Ex_delay;
E_Ex_delay = new unsigned int[E_Ex_Count];
for (unsigned int i=0;i<E_Ex_Count;++i)
E_Ex_delay[i]=vDelay.at(i);
for (int n=0;n<3;++n)
{
delete[] E_Ex_amp[n];

View File

@ -1,9 +1,98 @@
#include "engine.h"
#include "tools/AdrOp.h"
#include "tools/array_ops.h"
Engine::Engine()
Engine::Engine(Operator* op)
{
Op = op;
Init();
}
Engine::~Engine()
{
Reset();
}
void Engine::Init()
{
numTS = 0;
volt = Create_N_3DArray(Op->numLines);
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;iter<iterTS;++iter)
{
//voltage updates
for (pos[2]=1;pos[2]<Op->numLines[2]-1;++pos[2])
{
for (pos[1]=1;pos[1]<Op->numLines[1]-1;++pos[1])
{
for (pos[0]=1;pos[0]<Op->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;n<Op->E_Ex_Count;++n)
{
exc_pos = (int)numTS - (int)Op->E_Ex_delay[n];
exc_pos*= (exc_pos>0 && exc_pos<Op->ExciteLength);
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]<Op->numLines[0]-1;++pos[0])
{
for (pos[1]=0;pos[1]<Op->numLines[1]-1;++pos[1])
{
for (pos[2]=0;pos[2]<Op->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;
}

View File

@ -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

View File

@ -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<ExciteLength+1;++n)
{
ExciteSignal[n] = cos(2.0*PI*f0*(n*dT-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*n*dT/3.0-3,2));
// cerr << ExciteSignal[n] << endl;
}
}

View File

@ -19,10 +19,15 @@ public:
virtual int CalcECOperator() {};
virtual void CalcGaussianPulsExcitation(double f0, double fc);
virtual void ApplyElectricBC(bool* dirs) {}; //applied by default to all boundaries
virtual void ApplyMagneticBC(bool* dirs) {};
double GetTimestep() {return dT;};
double GetNumberCells();
void ShowSize();
virtual void Reset();
@ -35,18 +40,22 @@ protected:
unsigned int numLines[3];
//EC operator
FDTD_FLOAT* vv[3]; //calc new voltage from old voltage
FDTD_FLOAT* vi[3]; //calc new voltage from old current
FDTD_FLOAT* ii[3]; //calc new current from old current
FDTD_FLOAT* iv[3]; //calc new current from old voltage
FDTD_FLOAT**** vv; //calc new voltage from old voltage
FDTD_FLOAT**** vi; //calc new voltage from old current
FDTD_FLOAT**** ii; //calc new current from old current
FDTD_FLOAT**** iv; //calc new current from old voltage
//Excitation time-signal
unsigned int ExciteLength;
FDTD_FLOAT* ExciteSignal;
//E-Field Excitation
//! Calc the electric field excitation.
virtual bool CalcEFieldExcitation() {};
unsigned int E_Ex_Count;
unsigned int* E_Ex_index;
unsigned int* E_Ex_index[3];
FDTD_FLOAT* E_Ex_amp[3]; //represented as edge-voltages!!
FDTD_FLOAT* E_Ex_delay;
unsigned int* E_Ex_delay;
//Calc timestep only internal use
virtual double CalcTimestep() {};

View File

@ -1,6 +1,8 @@
#include <iostream>
#include <time.h>
#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);

View File

@ -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

58
tools/array_ops.cpp Normal file
View File

@ -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]<numLines[0];++pos[0])
{
array[pos[0]] = new FDTD_FLOAT*[numLines[1]];
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
{
array[pos[0]][pos[1]] = new FDTD_FLOAT[numLines[2]];
for (pos[2]=0;pos[2]<numLines[2];++pos[2])
{
array[pos[0]][pos[1]][pos[2]] = 0;
}
}
}
return array;
}
void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines)
{
if (array==NULL) return;
unsigned int pos[3];
for (pos[0]=0;pos[0]<numLines[0];++pos[0])
{
for (pos[1]=0;pos[1]<numLines[1];++pos[1])
{
delete[] array[pos[0]][pos[1]];
}
delete[] array[pos[0]];
}
delete[] array;
}
FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines)
{
FDTD_FLOAT**** array;
array = new FDTD_FLOAT***[3];
for (int n=0;n<3;++n)
{
array[n]=Create3DArray(numLines);
}
return array;
}
void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines)
{
if (array==NULL) return;
unsigned int pos[3];
for (int n=0;n<3;++n)
{
Delete3DArray(array[n],numLines);
}
delete[] array;
}

12
tools/array_ops.h Normal file
View File

@ -0,0 +1,12 @@
#ifndef ARRAY_OPS_H
#define ARRAY_OPS_H
#include "../FDTD/operator.h"
FDTD_FLOAT*** Create3DArray(unsigned int* numLines);
void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines);
FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines);
void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines);
#endif // ARRAY_OPS_H

View File

@ -3,5 +3,6 @@
#define __EPS0__ 8.85418781762e-12
#define __MUE0__ 1.256637062e-6
#define PI 3.141592653589793238462643383279
#endif // CONSTANTS_H