From b75476cc041f677aab1a429aca0c2220de0aa25a Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Sun, 15 Aug 2010 14:11:42 +0200 Subject: [PATCH 01/13] apply clear inheritance for cylindrical coords operator --- FDTD/operator_cylinder.cpp | 24 ++++++++++++------------ FDTD/operator_cylinder.h | 5 +---- 2 files changed, 13 insertions(+), 16 deletions(-) diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index f000030..ada236b 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -30,31 +30,31 @@ Operator_Cylinder* Operator_Cylinder::New(unsigned int numThreads) return op; } -Operator_Cylinder::Operator_Cylinder() : __OP_CYLINDER_BASE_CLASS__() +Operator_Cylinder::Operator_Cylinder() : Operator_Multithread() { m_MeshType = ProcessFields::CYLINDRICAL_MESH; } Operator_Cylinder::~Operator_Cylinder() { - __OP_CYLINDER_BASE_CLASS__::Reset(); + Operator_Multithread::Reset(); } void Operator_Cylinder::Init() { CC_closedAlpha = false; CC_R0_included = false; - __OP_CYLINDER_BASE_CLASS__::Init(); + Operator_Multithread::Init(); } void Operator_Cylinder::Reset() { - __OP_CYLINDER_BASE_CLASS__::Reset(); + Operator_Multithread::Reset(); } void Operator_Cylinder::InitOperator() { - __OP_CYLINDER_BASE_CLASS__::InitOperator(); + Operator_Multithread::InitOperator(); } inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const @@ -76,7 +76,7 @@ string Operator_Cylinder::GetDirName(int ny) const double Operator_Cylinder::GetMeshDelta(int n, const int* pos, bool dualMesh) const { - double delta = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(n,pos,dualMesh); + double delta = Operator_Multithread::GetMeshDelta(n,pos,dualMesh); if (delta==0) return delta; if (n==1) { @@ -89,7 +89,7 @@ double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) c { if (ny==2) { - double da = __OP_CYLINDER_BASE_CLASS__::GetMeshDelta(1,pos,dualMesh)/gridDelta; + double da = Operator_Multithread::GetMeshDelta(1,pos,dualMesh)/gridDelta; double r1,r2; if (!dualMesh) { @@ -105,12 +105,12 @@ double Operator_Cylinder::GetNodeArea(int ny, const int pos[3], bool dualMesh) c return da * pow(r2,2); return da/2* (pow(r2,2) - pow(r1,2)); } - return __OP_CYLINDER_BASE_CLASS__::GetNodeArea(ny,pos,dualMesh); + return Operator_Multithread::GetNodeArea(ny,pos,dualMesh); } bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) { - if (__OP_CYLINDER_BASE_CLASS__::SetGeometryCSX(geo)==false) return false; + if (Operator_Multithread::SetGeometryCSX(geo)==false) return false; double minmaxA = fabs(discLines[1][numLines[1]-1]-discLines[1][0]); if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width @@ -171,7 +171,7 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs) } } } - __OP_CYLINDER_BASE_CLASS__::ApplyElectricBC(dirs); + Operator_Multithread::ApplyElectricBC(dirs); } void Operator_Cylinder::ApplyMagneticBC(bool* dirs) @@ -185,7 +185,7 @@ void Operator_Cylinder::ApplyMagneticBC(bool* dirs) { dirs[0]=0; //no PMC in r_min directions... } - __OP_CYLINDER_BASE_CLASS__::ApplyMagneticBC(dirs); + Operator_Multithread::ApplyMagneticBC(dirs); } bool Operator_Cylinder::Calc_ECPos(int n, const unsigned int* pos, double* inEC) const @@ -398,7 +398,7 @@ bool Operator_Cylinder::Calc_ECPos(int n, const unsigned int* pos, double* inEC) bool Operator_Cylinder::Calc_EffMatPos(int n, const unsigned int* pos, double* inMat) const { - __OP_CYLINDER_BASE_CLASS__::Calc_EffMatPos(n, pos, inMat); + Operator_Multithread::Calc_EffMatPos(n, pos, inMat); // H_rho is not defined at position r==0 if (CC_R0_included && (n==0) && (pos[0]==0)) diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index ddf17e2..8e7e588 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -18,9 +18,6 @@ #ifndef OPERATOR_CYLINDER_H #define OPERATOR_CYLINDER_H -//! define the base class for the cylindrical coordinate FDTD -#define __OP_CYLINDER_BASE_CLASS__ Operator_Multithread - #include "operator_multithread.h" //! This class creates an operator for a cylindrical FDTD. @@ -28,7 +25,7 @@ This class creates an operator for a cylindrical FDTD. No special engine is necessary, all special cases e.g. a closed alpha mesh or an included r=0 case is treated by an operator/engine extension \sa operator_ext_cylinder. */ -class Operator_Cylinder : public __OP_CYLINDER_BASE_CLASS__ +class Operator_Cylinder : public Operator_Multithread { public: static Operator_Cylinder* New(unsigned int numThreads = 0); From e081a9cf94b2ba43fcb083fdeb670c069224231c Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Tue, 17 Aug 2010 10:26:38 +0200 Subject: [PATCH 02/13] Engine API: GetVolt/GetCurr methodes spilt up to Get/Set --- FDTD/engine.cpp | 16 +++++++++++++-- FDTD/engine.h | 13 ++++++++---- FDTD/engine_ext_cylinder.cpp | 20 +++++++++---------- FDTD/engine_ext_dispersive.cpp | 36 +++++++++++++++++----------------- FDTD/engine_ext_mur_abc.cpp | 12 ++++++------ FDTD/engine_ext_pml_sf.cpp | 12 ++++++------ FDTD/engine_sse.h | 13 ++++++++---- 7 files changed, 72 insertions(+), 50 deletions(-) diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index 147a332..b93aaca 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -125,13 +125,19 @@ void Engine::UpdateVoltages(unsigned int startX, unsigned int numX) void Engine::ApplyVoltageExcite() { int exc_pos; + unsigned int ny; + unsigned int pos[3]; //soft voltage excitation here (E-field excite) for (unsigned int n=0;nExc->Volt_Count;++n) { exc_pos = (int)numTS - (int)Op->Exc->Volt_delay[n]; exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length); // if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl; - GetVolt(Op->Exc->Volt_dir[n],Op->Exc->Volt_index[0][n],Op->Exc->Volt_index[1][n],Op->Exc->Volt_index[2][n]) += Op->Exc->Volt_amp[n]*Op->Exc->Signal_volt[exc_pos]; + ny = Op->Exc->Volt_dir[n]; + pos[0]=Op->Exc->Volt_index[0][n]; + pos[1]=Op->Exc->Volt_index[1][n]; + pos[2]=Op->Exc->Volt_index[2][n]; + SetVolt(ny,pos, GetVolt(ny,pos) + Op->Exc->Volt_amp[n]*Op->Exc->Signal_volt[exc_pos]); } // write the first excitation into the file "et" @@ -172,13 +178,19 @@ void Engine::UpdateCurrents(unsigned int startX, unsigned int numX) void Engine::ApplyCurrentExcite() { int exc_pos; + unsigned int ny; + unsigned int pos[3]; //soft current excitation here (H-field excite) for (unsigned int n=0;nExc->Curr_Count;++n) { exc_pos = (int)numTS - (int)Op->Exc->Curr_delay[n]; exc_pos *= (exc_pos>0 && exc_pos<=(int)Op->Exc->Length); // if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl; - GetCurr(Op->Exc->Curr_dir[n],Op->Exc->Curr_index[0][n],Op->Exc->Curr_index[1][n],Op->Exc->Curr_index[2][n]) += Op->Exc->Curr_amp[n]*Op->Exc->Signal_curr[exc_pos]; + ny = Op->Exc->Curr_dir[n]; + pos[0]=Op->Exc->Curr_index[0][n]; + pos[1]=Op->Exc->Curr_index[1][n]; + pos[2]=Op->Exc->Curr_index[2][n]; + SetCurr(ny,pos, GetCurr(ny,pos) + Op->Exc->Curr_amp[n]*Op->Exc->Signal_curr[exc_pos]); } // write the first excitation into the file "ht" diff --git a/FDTD/engine.h b/FDTD/engine.h index fd7acec..27d7493 100644 --- a/FDTD/engine.h +++ b/FDTD/engine.h @@ -47,10 +47,15 @@ public: virtual unsigned int GetNumberOfTimesteps() {return numTS;}; //this access functions muss be overloaded by any new engine using a different storage model - inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return volt[n][x][y][z]; } - inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int pos[3] ) { return volt[n][pos[0]][pos[1]][pos[2]]; } - inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return curr[n][x][y][z]; } - inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int pos[3] ) { return curr[n][pos[0]][pos[1]][pos[2]]; } + inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return volt[n][x][y][z]; } + inline virtual FDTD_FLOAT GetVolt( unsigned int n, const unsigned int pos[3] ) const { return volt[n][pos[0]][pos[1]][pos[2]]; } + inline virtual FDTD_FLOAT GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return curr[n][x][y][z]; } + inline virtual FDTD_FLOAT GetCurr( unsigned int n, const unsigned int pos[3] ) const { return curr[n][pos[0]][pos[1]][pos[2]]; } + + inline virtual void SetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value) { volt[n][x][y][z]=value; } + inline virtual void SetVolt( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value ) { volt[n][pos[0]][pos[1]][pos[2]]=value; } + inline virtual void SetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value) { curr[n][x][y][z]=value; } + inline virtual void SetCurr( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value ) { curr[n][pos[0]][pos[1]][pos[2]]=value; } virtual void UpdateVoltages(unsigned int startX, unsigned int numX); virtual void ApplyVoltageExcite(); diff --git a/FDTD/engine_ext_cylinder.cpp b/FDTD/engine_ext_cylinder.cpp index 1ac704e..3964791 100644 --- a/FDTD/engine_ext_cylinder.cpp +++ b/FDTD/engine_ext_cylinder.cpp @@ -40,18 +40,18 @@ void Engine_Ext_Cylinder::Apply2Voltages() pos[0] = 0; for (pos[2]=0;pos[2]GetVolt(2,0,0,pos[2]) *= cyl_Op->vv_R0[pos[2]]; + m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2])*cyl_Op->vv_R0[pos[2]]); for (pos[1]=0;pos[1]GetVolt(2,0,0,pos[2]) += cyl_Op->vi_R0[pos[2]] * m_Eng->GetCurr(1,0,pos[1],pos[2]); + m_Eng->SetVolt(2,0,0,pos[2], m_Eng->GetVolt(2,0,0,pos[2]) + cyl_Op->vi_R0[pos[2]] * m_Eng->GetCurr(1,0,pos[1],pos[2]) ); } } for (pos[1]=0;pos[1]GetVolt(1,0,pos[1],pos[2]) = 0; //no voltage in alpha-direction at r=0 - m_Eng->GetVolt(2,0,pos[1],pos[2]) = m_Eng->GetVolt(2,0,0,pos[2]); + m_Eng->SetVolt(1,0,pos[1],pos[2], 0); //no voltage in alpha-direction at r=0 + m_Eng->SetVolt(2,0,pos[1],pos[2], m_Eng->GetVolt(2,0,0,pos[2]) ); } } } @@ -64,9 +64,9 @@ void Engine_Ext_Cylinder::Apply2Voltages() { for (pos[2]=0;pos[2]GetVolt(0,pos[0],0,pos[2]) = m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]); - m_Eng->GetVolt(1,pos[0],0,pos[2]) = m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]); - m_Eng->GetVolt(2,pos[0],0,pos[2]) = m_Eng->GetVolt(2,pos[0],last_A_Line,pos[2]); + m_Eng->SetVolt(0,pos[0],0,pos[2], m_Eng->GetVolt(0,pos[0],last_A_Line,pos[2]) ); + m_Eng->SetVolt(1,pos[0],0,pos[2], m_Eng->GetVolt(1,pos[0],last_A_Line,pos[2]) ); + m_Eng->SetVolt(2,pos[0],0,pos[2], m_Eng->GetVolt(2,pos[0],last_A_Line,pos[2]) ); } } } @@ -83,9 +83,9 @@ void Engine_Ext_Cylinder::Apply2Current() unsigned int last_A_Line = numLines[1]-1; for (pos[2]=0;pos[2]GetCurr(0,pos[0],last_A_Line,pos[2]) = m_Eng->GetCurr(0,pos[0],0,pos[2]); - m_Eng->GetCurr(1,pos[0],last_A_Line,pos[2]) = m_Eng->GetCurr(1,pos[0],0,pos[2]); - m_Eng->GetCurr(2,pos[0],last_A_Line,pos[2]) = m_Eng->GetCurr(2,pos[0],0,pos[2]); + m_Eng->SetCurr(0,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(0,pos[0],0,pos[2]) ); + m_Eng->SetCurr(1,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(1,pos[0],0,pos[2]) ); + m_Eng->SetCurr(2,pos[0],last_A_Line,pos[2], m_Eng->GetCurr(2,pos[0],0,pos[2]) ); } } } diff --git a/FDTD/engine_ext_dispersive.cpp b/FDTD/engine_ext_dispersive.cpp index 8ee34fc..c2d4eb3 100644 --- a/FDTD/engine_ext_dispersive.cpp +++ b/FDTD/engine_ext_dispersive.cpp @@ -68,9 +68,9 @@ void Engine_Ext_Dispersive::Apply2Voltages() { for (unsigned int i=0;im_LM_Count;++i) { - m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[0][i]; - m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[1][i]; - m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[2][i]; + m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]); + m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]); + m_Eng->Engine::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]); } break; } @@ -79,18 +79,18 @@ void Engine_Ext_Dispersive::Apply2Voltages() Engine_sse* eng_sse = (Engine_sse*)m_Eng; for (unsigned int i=0;im_LM_Count;++i) { - eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[0][i]; - eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[1][i]; - eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[2][i]; + eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]); + eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]); + eng_sse->Engine_sse::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]); } break; } default: for (unsigned int i=0;im_LM_Count;++i) { - m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[0][i]; - m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[1][i]; - m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) -= volt_ADE[2][i]; + m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]); + m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]); + m_Eng->SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]); } break; } @@ -109,9 +109,9 @@ void Engine_Ext_Dispersive::Apply2Current() { for (unsigned int i=0;im_LM_Count;++i) { - m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[0][i]; - m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[1][i]; - m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[2][i]; + m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]); + m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]); + m_Eng->Engine::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]); } break; } @@ -120,18 +120,18 @@ void Engine_Ext_Dispersive::Apply2Current() Engine_sse* eng_sse = (Engine_sse*)m_Eng; for (unsigned int i=0;im_LM_Count;++i) { - eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[0][i]; - eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[1][i]; - eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[2][i]; + eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]); + eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]); + eng_sse->Engine_sse::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]); } break; } default: for (unsigned int i=0;im_LM_Count;++i) { - m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[0][i]; - m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[1][i]; - m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) -= curr_ADE[2][i]; + m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]); + m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]); + m_Eng->SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]); } break; } diff --git a/FDTD/engine_ext_mur_abc.cpp b/FDTD/engine_ext_mur_abc.cpp index b9370b1..f9da405 100644 --- a/FDTD/engine_ext_mur_abc.cpp +++ b/FDTD/engine_ext_mur_abc.cpp @@ -194,8 +194,8 @@ void Engine_Ext_Mur_ABC::Apply2Voltages() { for (pos[m_nyPP]=0;pos[m_nyPP]Engine::GetVolt(m_nyP,pos) = m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]; - m_Eng->Engine::GetVolt(m_nyPP,pos) = m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]; + m_Eng->Engine::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); + m_Eng->Engine::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); } } break; @@ -208,8 +208,8 @@ void Engine_Ext_Mur_ABC::Apply2Voltages() { for (pos[m_nyPP]=0;pos[m_nyPP]Engine_sse::GetVolt(m_nyP,pos) = m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]; - eng_sse->Engine_sse::GetVolt(m_nyPP,pos) = m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]; + eng_sse->Engine_sse::SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); + eng_sse->Engine_sse::SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); } } break; @@ -220,8 +220,8 @@ void Engine_Ext_Mur_ABC::Apply2Voltages() { for (pos[m_nyPP]=0;pos[m_nyPP]GetVolt(m_nyP,pos) = m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]; - m_Eng->GetVolt(m_nyPP,pos) = m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]; + m_Eng->SetVolt(m_nyP,pos, m_volt_nyP[pos[m_nyP]][pos[m_nyPP]]); + m_Eng->SetVolt(m_nyPP,pos, m_volt_nyPP[pos[m_nyP]][pos[m_nyPP]]); } } break; diff --git a/FDTD/engine_ext_pml_sf.cpp b/FDTD/engine_ext_pml_sf.cpp index 20beddd..d0e0271 100644 --- a/FDTD/engine_ext_pml_sf.cpp +++ b/FDTD/engine_ext_pml_sf.cpp @@ -181,9 +181,9 @@ void Engine_Ext_PML_SF_Plane::Apply2Voltages() for (pos[m_nyPP]=0;pos[m_nyPP]m_numLines[m_nyPP];++pos[m_nyPP]) { pml_pos[m_nyPP] = pos[m_nyPP]; - m_Eng->GetVolt(0,pos) = volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]]; - m_Eng->GetVolt(1,pos) = volt[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]]; - m_Eng->GetVolt(2,pos) = volt[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]]; + m_Eng->SetVolt(0,pos, volt[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); + m_Eng->SetVolt(1,pos, volt[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); + m_Eng->SetVolt(2,pos, volt[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] + volt[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); } } } @@ -254,9 +254,9 @@ void Engine_Ext_PML_SF_Plane::Apply2Current() { pml_pos[m_nyPP] = pos[m_nyPP]; - m_Eng->GetCurr(0,pos) = curr[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]]; - m_Eng->GetCurr(1,pos) = curr[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]]; - m_Eng->GetCurr(2,pos) = curr[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]]; + m_Eng->SetCurr(0,pos, curr[0][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][0][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); + m_Eng->SetCurr(1,pos, curr[0][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][1][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); + m_Eng->SetCurr(2,pos, curr[0][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] + curr[1][2][pml_pos[0]][pml_pos[1]][pml_pos[2]] ); } } } diff --git a/FDTD/engine_sse.h b/FDTD/engine_sse.h index 00cd324..7a53717 100644 --- a/FDTD/engine_sse.h +++ b/FDTD/engine_sse.h @@ -33,10 +33,15 @@ public: virtual unsigned int GetNumberOfTimesteps() {return numTS;}; //this access functions muss be overloaded by any new engine using a different storage model - inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return f4_volt[n][x][y][z%numVectors].f[z/numVectors]; } - inline virtual FDTD_FLOAT& GetVolt( unsigned int n, unsigned int pos[3] ) { return f4_volt[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; } - inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { return f4_curr[n][x][y][z%numVectors].f[z/numVectors]; } - inline virtual FDTD_FLOAT& GetCurr( unsigned int n, unsigned int pos[3] ) { return f4_curr[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; } + inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_volt[n][x][y][z%numVectors].f[z/numVectors]; } + inline virtual FDTD_FLOAT GetVolt( unsigned int n, const unsigned int pos[3] ) const { return f4_volt[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; } + inline virtual FDTD_FLOAT GetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_curr[n][x][y][z%numVectors].f[z/numVectors]; } + inline virtual FDTD_FLOAT GetCurr( unsigned int n, const unsigned int pos[3] ) const { return f4_curr[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]; } + + inline virtual void SetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value) { f4_volt[n][x][y][z%numVectors].f[z/numVectors]=value; } + inline virtual void SetVolt( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value ) { f4_volt[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]=value; } + inline virtual void SetCurr( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value) { f4_curr[n][x][y][z%numVectors].f[z/numVectors]=value; } + inline virtual void SetCurr( unsigned int n, const unsigned int pos[3], FDTD_FLOAT value ) { f4_curr[n][pos[0]][pos[1]][pos[2]%numVectors].f[pos[2]/numVectors]=value; } protected: Engine_sse(const Operator_sse* op); From 594b38e34560ce7463aee0b821304dfb2bbd7967 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 19 Aug 2010 10:12:09 +0200 Subject: [PATCH 03/13] operator: set a forced timestep --- FDTD/operator.cpp | 13 ++++++++++++- FDTD/operator.h | 2 ++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 9ca0926..77d333d 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -35,6 +35,7 @@ Operator::Operator() { m_MeshType = ProcessFields::CARTESIAN_MESH; Exc = 0; + dT = 0; } Operator::~Operator() @@ -77,6 +78,7 @@ void Operator::Init() m_BC[n]=0; Exc = 0; + dT = 0; } void Operator::Reset() @@ -531,7 +533,16 @@ int Operator::CalcECOperator() if (Calc_EC()==0) return -1; - CalcTimestep(); + if (dT>0) + { + double save_dT = dT; + CalcTimestep(); + if (dT Date: Thu, 19 Aug 2010 10:13:03 +0200 Subject: [PATCH 04/13] operator: setup excitation from xml --- FDTD/operator.h | 3 +++ openems.cpp | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/FDTD/operator.h b/FDTD/operator.h index fd3e43a..aac6788 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -25,6 +25,7 @@ class Operator_Extension; class Engine; +class TiXmlElement; //! Abstract base-class for the FDTD-operator class Operator @@ -44,6 +45,8 @@ public: virtual int CalcECOperator(); + virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return Exc->setupExcitation(Excite,maxTS);}; + inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; } inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; } diff --git a/openems.cpp b/openems.cpp index 23e8366..aa41423 100644 --- a/openems.cpp +++ b/openems.cpp @@ -351,7 +351,7 @@ int openEMS::SetupFDTD(const char* file) if ((maxTime_TS>0) && (maxTime_TSExc->setupExcitation( FDTD_Opts->FirstChildElement("Excitation"), NrTS )) + if (!FDTD_Op->SetupExcitation( FDTD_Opts->FirstChildElement("Excitation"), NrTS )) exit(2); if (DebugMat) From cabdf4a84ab8cc15bbbfae018508653030233684 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 2 Sep 2010 22:04:35 +0200 Subject: [PATCH 05/13] clear extensions method allows saver reset of MT engine --- FDTD/engine.cpp | 11 ++++++++--- FDTD/engine.h | 1 + FDTD/engine_multithread.cpp | 5 +++-- 3 files changed, 12 insertions(+), 5 deletions(-) diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index b93aaca..0d1e511 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -73,6 +73,13 @@ void Engine::InitExtensions() } } +void Engine::ClearExtensions() +{ + for (size_t n=0;n m_Eng_exts; ofstream file_et, file_ht; //excite signal dump file diff --git a/FDTD/engine_multithread.cpp b/FDTD/engine_multithread.cpp index b9c4c1d..55eed90 100644 --- a/FDTD/engine_multithread.cpp +++ b/FDTD/engine_multithread.cpp @@ -133,8 +133,9 @@ void Engine_Multithread::Init() void Engine_Multithread::Reset() { - if (!m_stopThreads) { - // prevent multiple invocations + if (!m_stopThreads) // prevent multiple invocations + { + ClearExtensions(); //prevent extensions from interfering with thread reset... // stop the threads //NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl; From 12c26f834a26bf7215e8b9ed76fe1862b4b0ed7c Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 2 Sep 2010 22:12:03 +0200 Subject: [PATCH 06/13] always show the optimal timestep --- FDTD/operator.cpp | 7 ++++++- FDTD/operator.h | 1 + 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 77d333d..1e87285 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -304,7 +304,10 @@ void Operator::ShowStat() const cout << "in " << GetDirName(1) << " direction\t\t: " << m_Nr_PEC[1] << endl; cout << "in " << GetDirName(2) << " direction\t\t: " << m_Nr_PEC[2] << endl; cout << "-----------------------------------" << endl; - cout << "Timestep (s)\t\t: " << dT << endl; + cout << "Timestep (s)\t\t: " << dT ; + if (opt_dT) + cout <<"\t(" << opt_dT << ")"; + cout << endl; cout << "Timestep method name\t: " << m_Used_TS_Name << endl; cout << "Nyquist criteria (TS)\t: " << Exc->GetNyquistNum() << endl; cout << "Nyquist criteria (s)\t: " << Exc->GetNyquistNum()*dT << endl; @@ -533,10 +536,12 @@ int Operator::CalcECOperator() if (Calc_EC()==0) return -1; + opt_dT = 0; if (dT>0) { double save_dT = dT; CalcTimestep(); + opt_dT = dT; if (dT Date: Thu, 2 Sep 2010 22:14:40 +0200 Subject: [PATCH 07/13] missing gpl header in cylinder extension --- FDTD/operator_ext_cylinder.cpp | 17 +++++++++++++++++ FDTD/operator_ext_cylinder.h | 17 +++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/FDTD/operator_ext_cylinder.cpp b/FDTD/operator_ext_cylinder.cpp index e29eae1..4c3fbf9 100644 --- a/FDTD/operator_ext_cylinder.cpp +++ b/FDTD/operator_ext_cylinder.cpp @@ -1,3 +1,20 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + #include "operator_ext_cylinder.h" #include "operator_cylinder.h" #include "engine_ext_cylinder.h" diff --git a/FDTD/operator_ext_cylinder.h b/FDTD/operator_ext_cylinder.h index 26730d2..b880d39 100644 --- a/FDTD/operator_ext_cylinder.h +++ b/FDTD/operator_ext_cylinder.h @@ -1,3 +1,20 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + #ifndef OPERATOR_EXT_CYLINDER_H #define OPERATOR_EXT_CYLINDER_H From db0f4ab3e09a174eda80e6763515acebb95ac877 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 2 Sep 2010 22:16:37 +0200 Subject: [PATCH 08/13] MT-engine has access to MT-operator --- FDTD/engine_multithread.cpp | 3 ++- FDTD/engine_multithread.h | 3 ++- FDTD/operator_multithread.h | 1 + 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/FDTD/engine_multithread.cpp b/FDTD/engine_multithread.cpp index 55eed90..44b4121 100644 --- a/FDTD/engine_multithread.cpp +++ b/FDTD/engine_multithread.cpp @@ -44,8 +44,9 @@ Engine_Multithread* Engine_Multithread::New(const Operator_Multithread* op, unsi return e; } -Engine_Multithread::Engine_Multithread(const Operator_SSE_Compressed* op) : Engine_SSE_Compressed(op) +Engine_Multithread::Engine_Multithread(const Operator_Multithread* op) : Engine_SSE_Compressed(op) { + m_Op_MT = op; m_type = SSE; m_barrier_VoltUpdate = 0; m_barrier_VoltExcite = 0; diff --git a/FDTD/engine_multithread.h b/FDTD/engine_multithread.h index 6b5ffe3..07f9ec6 100644 --- a/FDTD/engine_multithread.h +++ b/FDTD/engine_multithread.h @@ -90,7 +90,8 @@ public: virtual bool IterateTS(unsigned int iterTS); protected: - Engine_Multithread(const Operator_SSE_Compressed* op); + Engine_Multithread(const Operator_Multithread* op); + const Operator_Multithread* m_Op_MT; boost::thread_group m_thread_group; boost::barrier *m_startBarrier, *m_stopBarrier; boost::barrier *m_barrier_VoltUpdate, *m_barrier_VoltExcite, *m_barrier_PreVolt, *m_barrier_PostVolt; diff --git a/FDTD/operator_multithread.h b/FDTD/operator_multithread.h index ddc17af..d6333d2 100644 --- a/FDTD/operator_multithread.h +++ b/FDTD/operator_multithread.h @@ -24,6 +24,7 @@ class Operator_Multithread : public Operator_SSE_Compressed { + friend class Engine_Multithread; friend class Operator_Thread; public: //! Create a new operator From 1a818f659b8add0cf23a7ab023d5f18ac68d3d98 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Fri, 3 Sep 2010 11:36:59 +0200 Subject: [PATCH 09/13] Read a forced timestep from xml --- openems.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/openems.cpp b/openems.cpp index aa41423..bf69279 100644 --- a/openems.cpp +++ b/openems.cpp @@ -345,6 +345,11 @@ int openEMS::SetupFDTD(const char* file) if (CSX.GetQtyPropertyType(CSProperties::LORENTZMATERIAL)>0) FDTD_Op->AddExtension(new Operator_Ext_LorentzMaterial(FDTD_Op)); + double timestep=0; + FDTD_Opts->QueryDoubleAttribute("TimeStep",×tep); + if (timestep) + FDTD_Op->SetTimestep(timestep); + FDTD_Op->CalcECOperator(); unsigned int maxTime_TS = (unsigned int)(maxTime/FDTD_Op->GetTimestep()); From fc2b60ba3ed181cb000b6eb5ca34f94aca84c2e8 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Fri, 3 Sep 2010 11:53:33 +0200 Subject: [PATCH 10/13] Operator Extension clone method: allows to create a copy/clone of an existing extension This will be necessary for the upcoming multi-grid approach... --- FDTD/operator_ext_mur_abc.cpp | 36 +++++++++++++++++++++++++++-------- FDTD/operator_ext_mur_abc.h | 5 +++++ FDTD/operator_extension.cpp | 13 +++++++++++++ FDTD/operator_extension.h | 9 +++++++++ 4 files changed, 55 insertions(+), 8 deletions(-) diff --git a/FDTD/operator_ext_mur_abc.cpp b/FDTD/operator_ext_mur_abc.cpp index d1f37f6..327c587 100644 --- a/FDTD/operator_ext_mur_abc.cpp +++ b/FDTD/operator_ext_mur_abc.cpp @@ -21,6 +21,33 @@ #include "tools/array_ops.h" Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op) +{ + Initialize(); +} + +Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC() +{ + Delete2DArray(m_Mur_Coeff_nyP,m_numLines); + m_Mur_Coeff_nyP = NULL; + Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); + m_Mur_Coeff_nyPP = NULL; +} + +Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op, Operator_Ext_Mur_ABC* op_ext) : Operator_Extension(op, op_ext) +{ + Initialize(); + m_v_phase = op_ext->m_v_phase; + SetDirection(op_ext->m_ny,op_ext->m_top); +} + +Operator_Extension* Operator_Ext_Mur_ABC::Clone(Operator* op) +{ + if (dynamic_cast(this)==NULL) + return NULL; + return new Operator_Ext_Mur_ABC(op, this); +} + +void Operator_Ext_Mur_ABC::Initialize() { m_ny = -1; m_nyP = -1; @@ -37,14 +64,6 @@ Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op m_numLines[1]=0; } -Operator_Ext_Mur_ABC::~Operator_Ext_Mur_ABC() -{ - Delete2DArray(m_Mur_Coeff_nyP,m_numLines); - m_Mur_Coeff_nyP = NULL; - Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); - m_Mur_Coeff_nyPP = NULL; -} - void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny) { if ((ny<0) || (ny>2)) @@ -54,6 +73,7 @@ void Operator_Ext_Mur_ABC::SetDirection(int ny, bool top_ny) Delete2DArray(m_Mur_Coeff_nyPP,m_numLines); m_ny = ny; + m_top = top_ny; m_nyP = (ny+1)%3; m_nyPP = (ny+2)%3; if (!top_ny) diff --git a/FDTD/operator_ext_mur_abc.h b/FDTD/operator_ext_mur_abc.h index 0d6b340..798886c 100644 --- a/FDTD/operator_ext_mur_abc.h +++ b/FDTD/operator_ext_mur_abc.h @@ -28,6 +28,9 @@ public: Operator_Ext_Mur_ABC(Operator* op); ~Operator_Ext_Mur_ABC(); + Operator_Ext_Mur_ABC(Operator* op, Operator_Ext_Mur_ABC* op_ext); + virtual Operator_Extension* Clone(Operator* op); + //! Define the direction of this ABC: \a ny=0,1,2 -> x,y,z and if at bottom_ny -> e.g. x=0 or x=end void SetDirection(int ny, bool top_ny); @@ -45,8 +48,10 @@ public: virtual void ShowStat(ostream &ostr) const; protected: + void Initialize(); int m_ny; int m_nyP,m_nyPP; + bool m_top; unsigned int m_LineNr; int m_LineNr_Shift; diff --git a/FDTD/operator_extension.cpp b/FDTD/operator_extension.cpp index 69f0442..cc3b3bc 100644 --- a/FDTD/operator_extension.cpp +++ b/FDTD/operator_extension.cpp @@ -24,7 +24,20 @@ Operator_Extension::Operator_Extension(Operator* op) m_Op = op; } +Operator_Extension::Operator_Extension(Operator* op, Operator_Extension* op_ext) +{ + UNUSED(op_ext); + m_Op = op; +} + void Operator_Extension::ShowStat(ostream &ostr) const { ostr << "--- " << GetExtensionName() << " ---" << endl; } + +Operator_Extension* Operator_Extension::Clone(Operator* op) +{ + if (dynamic_cast(this)==NULL) + return NULL; + return new Operator_Extension(op, this); +} diff --git a/FDTD/operator_extension.h b/FDTD/operator_extension.h index fcc24cb..54d24d9 100644 --- a/FDTD/operator_extension.h +++ b/FDTD/operator_extension.h @@ -35,6 +35,13 @@ class Operator_Extension { friend class Engine_Extension; public: + //! Create a clone of this extension, will return NULL if this is impossible + /*! + Create a clone of this extension, will return NULL if this is impossible (e.g. derived extension has no clone method and copy-constructor)... + BuildExtension has to be called separatly! + */ + virtual Operator_Extension* Clone(Operator* op); + virtual bool BuildExtension() {return true;} virtual Engine_Extension* CreateEngineExtention() {return 0;} @@ -47,6 +54,8 @@ public: protected: Operator_Extension(Operator* op); + //! Copy constructor, returns NULL if extension cannot be copied... + Operator_Extension(Operator* op, Operator_Extension* op_ext); Operator* m_Op; }; From a52cd4711af57f9bf8ab7b7cbecd4172b71c9b36 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Fri, 3 Sep 2010 12:14:25 +0200 Subject: [PATCH 11/13] MT operator: new separate calc start-stop lines method --- FDTD/engine_multithread.cpp | 13 ++++++++----- FDTD/operator_multithread.cpp | 31 ++++++++++++++++++++++--------- FDTD/operator_multithread.h | 8 ++++++++ 3 files changed, 38 insertions(+), 14 deletions(-) diff --git a/FDTD/engine_multithread.cpp b/FDTD/engine_multithread.cpp index 44b4121..2e78bc2 100644 --- a/FDTD/engine_multithread.cpp +++ b/FDTD/engine_multithread.cpp @@ -96,6 +96,10 @@ void Engine_Multithread::Init() if (m_numThreads == 0) m_numThreads = boost::thread::hardware_concurrency(); + vector m_Start_Lines; + vector m_Stop_Lines; + m_Op_MT->CalcStartStopLines( m_numThreads, m_Start_Lines, m_Stop_Lines ); + cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: ("; m_barrier_VoltUpdate = new boost::barrier(m_numThreads); // numThread workers m_barrier_VoltExcite = new boost::barrier(m_numThreads+1); // numThread workers + 1 excitation thread @@ -110,15 +114,14 @@ void Engine_Multithread::Init() m_startBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller m_stopBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller - unsigned int linesPerThread = round((float)numLines[0] / (float)m_numThreads); for (unsigned int n=0; n &start, vector &stop) const +{ + unsigned int linesPerThread = round((float)numLines[0] / (float)numThreads); + if ((numThreads-1) * linesPerThread >= numLines[0]) + --numThreads; + + start.resize(numThreads); + stop.resize(numThreads); + + for (unsigned int n=0; n= numLines[0]) - --m_numThreads; + vector m_Start_Lines; + vector m_Stop_Lines; + CalcStartStopLines( m_numThreads, m_Start_Lines, m_Stop_Lines ); cout << "Multithreaded operator using " << m_numThreads << " threads." << std::endl; @@ -94,12 +112,7 @@ int Operator_Multithread::CalcECOperator() for (unsigned int n=0; n &start, vector &stop) const; }; class Operator_Thread From bd4794ecc4d0438cb6a5cf4ffb24924357d808df Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Wed, 8 Sep 2010 07:36:32 +0200 Subject: [PATCH 12/13] NEW: cylindrical multigrid FDTD this is a new multi grid approach for the cylindrical FDTD. The FDTD domain will be split in two regions in radial direction. The "inner" region will have half as many disc-lines in alpha direction and therefore allow for a much larger timestep which increases the simulation speed. Todo: - currently only a homogeneous disc is allowed in alpha direction - some extensions have to be tested and prepared for this approach (e.g. pml) - speed enhancement and more efficient memory usage - lots and lots of testing... --- FDTD/engine_cylindermultigrid.cpp | 250 ++++++++++++++++++++++++++ FDTD/engine_cylindermultigrid.h | 82 +++++++++ FDTD/engine_ext_cylindermultigrid.cpp | 158 ++++++++++++++++ FDTD/engine_ext_cylindermultigrid.h | 58 ++++++ FDTD/engine_multithread.h | 1 + FDTD/engine_sse.h | 1 + FDTD/operator.cpp | 8 +- FDTD/operator.h | 2 +- FDTD/operator_cylinder.h | 1 + FDTD/operator_cylindermultigrid.cpp | 189 +++++++++++++++++++ FDTD/operator_cylindermultigrid.h | 70 ++++++++ FDTD/operator_multithread.cpp | 3 + openEMS.pro | 10 +- openems.cpp | 12 +- 14 files changed, 836 insertions(+), 9 deletions(-) create mode 100644 FDTD/engine_cylindermultigrid.cpp create mode 100644 FDTD/engine_cylindermultigrid.h create mode 100644 FDTD/engine_ext_cylindermultigrid.cpp create mode 100644 FDTD/engine_ext_cylindermultigrid.h create mode 100644 FDTD/operator_cylindermultigrid.cpp create mode 100644 FDTD/operator_cylindermultigrid.h diff --git a/FDTD/engine_cylindermultigrid.cpp b/FDTD/engine_cylindermultigrid.cpp new file mode 100644 index 0000000..c517f4e --- /dev/null +++ b/FDTD/engine_cylindermultigrid.cpp @@ -0,0 +1,250 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#include "engine_cylindermultigrid.h" +#include "operator_cylindermultigrid.h" +#include "engine_ext_cylindermultigrid.h" + +Engine_CylinderMultiGrid* Engine_CylinderMultiGrid::New(const Operator_CylinderMultiGrid* op, unsigned int numThreads) +{ + cout << "Create FDTD engine (cylindrical multi grid mesh using sse compression + multithreading)" << endl; + Engine_CylinderMultiGrid* e = new Engine_CylinderMultiGrid(op); + e->setNumThreads( numThreads ); + e->Init(); + return e; +} + +Engine_CylinderMultiGrid::Engine_CylinderMultiGrid(const Operator_CylinderMultiGrid* op) : Engine_Multithread(op) +{ + Op_CMG = op; + + m_WaitOnBase = new boost::barrier(2); + m_WaitOnChild = new boost::barrier(2); + m_WaitOnSync = new boost::barrier(2); + + Engine* eng = op->GetInnerOperator()->CreateEngine(); + m_InnerEngine = dynamic_cast(eng); + + m_Eng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,true); + m_Eng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync); + m_Eng_Ext_MG->SetEngine(this); + m_InnerEng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,false); + m_InnerEng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync); + + m_InnerEngine->m_Eng_exts.push_back(m_InnerEng_Ext_MG); +} + +Engine_CylinderMultiGrid::~Engine_CylinderMultiGrid() +{ + m_Thread_NumTS = 0; + m_startBarrier->wait(); + + m_IteratorThread_Group.join_all(); + + delete m_InnerEngine; + m_InnerEngine = NULL; + + delete m_WaitOnBase; + m_WaitOnBase = NULL; + delete m_WaitOnChild; + m_WaitOnChild = NULL; + delete m_WaitOnSync; + m_WaitOnSync = NULL; + + delete m_startBarrier; + m_startBarrier = NULL; + delete m_stopBarrier; + m_stopBarrier = NULL; + +} + +void Engine_CylinderMultiGrid::Init() +{ + Engine_Multithread::Init(); + + m_Eng_exts.push_back(m_Eng_Ext_MG); + + m_startBarrier = new boost::barrier(3); //both engines + organizer + m_stopBarrier = new boost::barrier(3); //both engines + organizer + + boost::thread *t = NULL; + + t = new boost::thread( Engine_CylinderMultiGrid_Thread(this,m_startBarrier,m_stopBarrier,&m_Thread_NumTS) ); + m_IteratorThread_Group.add_thread( t ); + + t = new boost::thread( Engine_CylinderMultiGrid_Thread(m_InnerEngine,m_startBarrier,m_stopBarrier,&m_Thread_NumTS) ); + m_IteratorThread_Group.add_thread( t ); +} + +bool Engine_CylinderMultiGrid::IterateTS(unsigned int iterTS) +{ + m_Thread_NumTS = iterTS; + + m_startBarrier->wait(); //start base and child iterations + + m_stopBarrier->wait(); //tell base and child to wait for another start event... + + //interpolate child data to base mesh... + for (unsigned int n=0;nm_Split_Pos-1;++n) + InterpolVoltChild2Base(n); + for (unsigned int n=0;nm_Split_Pos-2;++n) + InterpolCurrChild2Base(n); + + return true; +} + +void Engine_CylinderMultiGrid::InitExtensions() +{ + m_InnerEngine->InitExtensions(); + Engine_Multithread::InitExtensions(); +} + +void Engine_CylinderMultiGrid::InterpolVoltChild2Base(unsigned int rzPlane) +{ + //interpolate voltages from child engine to the base engine... + unsigned int pos[3]; + unsigned int pos1_rz1, pos1_rz2; + unsigned int pos1_a1, pos1_a2; + pos[0] = rzPlane; + bool isOdd, isEven; + f4vector half, one_eighth, three_eighth; + for (int n=0;n<4;++n) + { + half.f[n]=0.5; + one_eighth.f[n] = 1.0/8.0; + three_eighth.f[n] = 3.0/8.0; + } + for (pos[1]=0; pos[1]0); + + pos1_rz1 = pos[1]/2; + pos1_rz2 = pos[1]/2 + isOdd; + + pos1_rz1 += pos[1]==1; + pos1_rz2 -= pos[1]==(numLines[1]-2); + + pos1_a1 = pos[1]/2; + pos1_a2 = pos[1]/2 + isOdd - isEven; + pos1_a2 += pos[1]==(numLines[1]-1); + pos1_a2 -= pos[1]==(numLines[1]-2); + + for (pos[2]=0; pos[2]f4_volt[0][pos[0]][pos1_rz1][pos[2]].v + m_InnerEngine->f4_volt[0][pos[0]][pos1_rz2][pos[2]].v); + + //z - direction + f4_volt[2][pos[0]][pos[1]][pos[2]].v = half.v * ( m_InnerEngine->f4_volt[2][pos[0]][pos1_rz1][pos[2]].v + m_InnerEngine->f4_volt[2][pos[0]][pos1_rz2][pos[2]].v); + + //alpha - direction + f4_volt[1][pos[0]][pos[1]][pos[2]].v = three_eighth.v * m_InnerEngine->f4_volt[1][pos[0]][pos1_a1][pos[2]].v; + f4_volt[1][pos[0]][pos[1]][pos[2]].v += one_eighth.v * m_InnerEngine->f4_volt[1][pos[0]][pos1_a2][pos[2]].v; + } + } + + //r,z - interpolation correction... + pos[1]=1; + for (pos[2]=0; pos[2]f4_volt[0][pos[0]][1][pos[2]].v - m_InnerEngine->f4_volt[0][pos[0]][2][pos[2]].v); + f4_volt[2][pos[0]][pos[1]][pos[2]].v += half.v * (m_InnerEngine->f4_volt[2][pos[0]][1][pos[2]].v - m_InnerEngine->f4_volt[2][pos[0]][2][pos[2]].v); + } + pos[1]=numLines[1]-2; + for (pos[2]=0; pos[2]f4_volt[0][pos[0]][pos[1]/2][pos[2]].v - m_InnerEngine->f4_volt[0][pos[0]][pos[1]/2-1][pos[2]].v); + f4_volt[2][pos[0]][pos[1]][pos[2]].v += half.v * (m_InnerEngine->f4_volt[2][pos[0]][pos[1]/2][pos[2]].v - m_InnerEngine->f4_volt[2][pos[0]][pos[1]/2-1][pos[2]].v); + } +} + +void Engine_CylinderMultiGrid::InterpolCurrChild2Base(unsigned int rzPlane) +{ + //interpolate voltages from child engine to the base engine... + unsigned int pos[3]; + unsigned int pos1_rz1, pos1_rz2; + unsigned int pos1_a1, pos1_a2; + pos[0] = rzPlane; + bool isOdd, isEven; + f4vector quarter, one_fourth, three_fourth; + for (int n=0;n<4;++n) + { + quarter.f[n]=0.25; + one_fourth.f[n] = 1.0/4.0; + three_fourth.f[n] = 3.0/4.0; + } + for (pos[1]=0; pos[1]0); + + pos1_a1 = pos[1]/2; + pos1_a2 = pos[1]/2 + isOdd; + + pos1_a1 += pos[1]==1; + pos1_a2 -= pos[1]==(numLines[1]-2); + + pos1_rz1 = pos[1]/2; + pos1_rz2 = pos[1]/2 + isOdd - isEven; + pos1_rz2 += pos[1]==(numLines[1]-1); + pos1_rz2 -= pos[1]==(numLines[1]-2); + + for (pos[2]=0; pos[2]f4_curr[0][pos[0]][pos1_rz1][pos[2]].v; + f4_curr[0][pos[0]][pos[1]][pos[2]].v += one_fourth.v * m_InnerEngine->f4_curr[0][pos[0]][pos1_rz2][pos[2]].v; + + //z - direction + f4_curr[2][pos[0]][pos[1]][pos[2]].v = three_fourth.v * m_InnerEngine->f4_curr[2][pos[0]][pos1_rz1][pos[2]].v; + f4_curr[2][pos[0]][pos[1]][pos[2]].v += one_fourth.v * m_InnerEngine->f4_curr[2][pos[0]][pos1_rz2][pos[2]].v; + + //alpha - direction + f4_curr[1][pos[0]][pos[1]][pos[2]].v = quarter.v * (m_InnerEngine->f4_curr[1][pos[0]][pos1_a1][pos[2]].v + m_InnerEngine->f4_curr[1][pos[0]][pos1_a2][pos[2]].v); + } + + } + //alpha - interpolation correction... + pos[1]=1; + for (pos[2]=0; pos[2]f4_curr[1][pos[0]][1][pos[2]].v - m_InnerEngine->f4_curr[1][pos[0]][2][pos[2]].v); + pos[1]=numLines[1]-2; + for (pos[2]=0; pos[2]f4_curr[1][pos[0]][pos[1]/2][pos[2]].v - m_InnerEngine->f4_curr[1][pos[0]][pos[1]/2-1][pos[2]].v); +} + +/****************************************************************************************/ +Engine_CylinderMultiGrid_Thread::Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS) +{ + m_startBarrier = start; + m_stopBarrier = stop; + m_Eng=engine; + m_numTS = numTS; +} + +void Engine_CylinderMultiGrid_Thread::operator()() +{ + m_startBarrier->wait(); //wait for Base engine to start the iterations... + + while(*m_numTS>0) //m_numTS==0 request to terminate this thread... + { + m_Eng->Engine_Multithread::IterateTS(*m_numTS); + m_stopBarrier->wait(); //sync all workers after iterations are performed + m_startBarrier->wait(); //wait for Base engine to start the iterations again ... + } +} diff --git a/FDTD/engine_cylindermultigrid.h b/FDTD/engine_cylindermultigrid.h new file mode 100644 index 0000000..366f7b4 --- /dev/null +++ b/FDTD/engine_cylindermultigrid.h @@ -0,0 +1,82 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#ifndef ENGINE_CYLINDERMULTIGRID_H +#define ENGINE_CYLINDERMULTIGRID_H + +#include "engine_multithread.h" + +class Operator_CylinderMultiGrid; +class Engine_CylinderMultiGrid_Thread; +class Engine_Ext_CylinderMultiGrid; + +class Engine_CylinderMultiGrid : public Engine_Multithread +{ + friend class Engine_Ext_CylinderMultiGrid; +public: + Engine_CylinderMultiGrid(); + + static Engine_CylinderMultiGrid* New(const Operator_CylinderMultiGrid* op, unsigned int numThreads = 0); + virtual ~Engine_CylinderMultiGrid(); + + virtual void InterpolVoltChild2Base(unsigned int rzPlane); + virtual void InterpolCurrChild2Base(unsigned int rzPlane); + + virtual void Init(); + + //! Iterate \a iterTS number of timesteps + virtual bool IterateTS(unsigned int iterTS); + + virtual void InitExtensions(); + +protected: + Engine_CylinderMultiGrid(const Operator_CylinderMultiGrid* op); + const Operator_CylinderMultiGrid* Op_CMG; + + Engine_Multithread* m_InnerEngine; + + volatile unsigned int m_Thread_NumTS; + boost::thread_group m_IteratorThread_Group; + boost::barrier *m_startBarrier; + boost::barrier *m_stopBarrier; + Engine_CylinderMultiGrid_Thread* m_IteratorThread; + Engine_CylinderMultiGrid_Thread* m_InnerIteratorThread; + + //extension barrier + boost::barrier *m_WaitOnBase; + boost::barrier *m_WaitOnChild; + boost::barrier *m_WaitOnSync; + + Engine_Ext_CylinderMultiGrid* m_Eng_Ext_MG; + Engine_Ext_CylinderMultiGrid* m_InnerEng_Ext_MG; +}; + + +class Engine_CylinderMultiGrid_Thread +{ +public: + Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS); + void operator()(); + +protected: + Engine_Multithread *m_Eng; + boost::barrier *m_startBarrier; + boost::barrier *m_stopBarrier; + volatile unsigned int *m_numTS; +}; + +#endif // ENGINE_CYLINDERMULTIGRID_H diff --git a/FDTD/engine_ext_cylindermultigrid.cpp b/FDTD/engine_ext_cylindermultigrid.cpp new file mode 100644 index 0000000..e8d858d --- /dev/null +++ b/FDTD/engine_ext_cylindermultigrid.cpp @@ -0,0 +1,158 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#include "engine.h" +#include "engine_ext_cylindermultigrid.h" +#include "engine_cylindermultigrid.h" + +Engine_Ext_CylinderMultiGrid::Engine_Ext_CylinderMultiGrid(Operator_Extension* op_ext, bool isBase) : Engine_Extension(op_ext) +{ + m_IsBase = isBase; + m_Eng_MG = NULL; +} + +Engine_Ext_CylinderMultiGrid::~Engine_Ext_CylinderMultiGrid() +{ +} + +void Engine_Ext_CylinderMultiGrid::SetBarrier(boost::barrier* waitBase, boost::barrier* waitChild, boost::barrier* waitSync) +{ + m_WaitOnBase = waitBase; + m_WaitOnChild = waitChild; + m_WaitOnSync = waitSync; +} + +void Engine_Ext_CylinderMultiGrid::SetEngine(Engine* eng) +{ + m_Eng_MG = dynamic_cast(eng); + if (m_Eng_MG==NULL) + { + cerr << "Engine_Ext_CylinderMultiGrid::SetEngine(): Error" << endl; + exit(0); + } +} + +void Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates() +{ + //cerr << "Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates() for " << m_IsBase << endl; + if (!m_IsBase) + { + //cerr << "child: volt wait on base " << endl; + m_WaitOnBase->wait(); //wait on base to finisch current sync and/or to finisch voltage updates, than start child voltage updates + } +} + +void Engine_Ext_CylinderMultiGrid::DoPostVoltageUpdates() +{ + +} + +void Engine_Ext_CylinderMultiGrid::Apply2Voltages() +{ + if (m_IsBase) + { + m_WaitOnBase->wait(); //base voltage updates are done, tell child to start its voltage updates + m_WaitOnChild->wait(); //wait for child to finisch its updates + SyncVoltages(); //child is finisch, run sync and go to current updates next + m_WaitOnSync->wait(); //sync is done... move on and tell child to move on... + } + else + { + m_WaitOnChild->wait(); //child is finished voltage updates, will tell base to run sync + m_WaitOnSync->wait(); //wait for base to finisch sync before going to wait for current updates + } +} + +void Engine_Ext_CylinderMultiGrid::SyncVoltages() +{ + if (m_Eng_MG==NULL) + { + cerr << "Engine_Ext_CylinderMultiGrid::SyncVoltages: Error engine is NULL" << endl; + return; + } + + unsigned int* numLines = m_Eng_MG->numLines; + + Engine_Multithread* m_InnerEng = m_Eng_MG->m_InnerEngine; + + //interpolate voltages from base engine to child engine... + unsigned int pos[3]; + pos[0] = m_Eng_MG->Op_CMG->GetSplitPos()-1; + unsigned int pos1_half = 0; + f4vector v_null; v_null.f[0] = 0;v_null.f[1] = 0;v_null.f[2] = 0;v_null.f[3] = 0; + for (pos[1]=0; pos[1]numVectors; ++pos[2]) + { + //r - direczion + m_InnerEng->f4_volt[0][pos[0]][pos1_half][pos[2]].v = v_null.v; + + //z - direction + m_InnerEng->f4_volt[2][pos[0]][pos1_half][pos[2]].v = m_Eng_MG->f4_volt[2][pos[0]][pos[1]][pos[2]].v; + + //alpha - direction + m_InnerEng->f4_volt[1][pos[0]][pos1_half][pos[2]].v = m_Eng_MG->f4_volt[1][pos[0]][pos[1]][pos[2]].v; + m_InnerEng->f4_volt[1][pos[0]][pos1_half][pos[2]].v += m_Eng_MG->f4_volt[1][pos[0]][pos[1]+1][pos[2]].v; + } + } +} + +void Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates() +{ + //cerr << "Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates() for " << m_IsBase << endl; + if (!m_IsBase) + { + //cerr << "child: curr wait on base " << endl; + m_WaitOnBase->wait(); //wait on base to finisch voltage sync and current updates, than start child current updates + } +} + +void Engine_Ext_CylinderMultiGrid::DoPostCurrentUpdates() +{ + +} + +void Engine_Ext_CylinderMultiGrid::Apply2Current() +{ + if (m_IsBase) + { + //cerr << "Base: curr wait on base done, wait on sync" << endl; + m_WaitOnBase->wait(); //base current updates are done, tell child to start its current updates + m_WaitOnChild->wait(); //wait for child to finisch its updates + SyncCurrents(); //child is finisch, run sync and go to voltage updates next + m_WaitOnSync->wait(); //sync is done... move on and tell child to move on... + } + else + { + m_WaitOnChild->wait(); //child is finished current updates, will tell base to run sync... + m_WaitOnSync->wait(); //wait for base to finisch sync before going to wait for next voltage updates + //cerr << "Child: curr done, wait on sync" << endl; + } +} + +void Engine_Ext_CylinderMultiGrid::SyncCurrents() +{ + if (m_Eng_MG==NULL) + { + cerr << "Engine_Ext_CylinderMultiGrid::SyncCurrents: Error engine is NULL" << endl; + return; + } + + m_Eng_MG->InterpolCurrChild2Base(m_Eng_MG->Op_CMG->GetSplitPos()-2); + return; +} diff --git a/FDTD/engine_ext_cylindermultigrid.h b/FDTD/engine_ext_cylindermultigrid.h new file mode 100644 index 0000000..b995a1c --- /dev/null +++ b/FDTD/engine_ext_cylindermultigrid.h @@ -0,0 +1,58 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#ifndef ENGINE_EXT_CYLINDERMULTIGRID_H +#define ENGINE_EXT_CYLINDERMULTIGRID_H + +#include "engine_cylindermultigrid.h" +#include "engine_extension.h" +#include "operator_cylindermultigrid.h" + +class Operator_Ext_CylinderMultiGrid; + +class Engine_Ext_CylinderMultiGrid : public Engine_Extension +{ +public: + Engine_Ext_CylinderMultiGrid(Operator_Extension* op_ext, bool isBase); + virtual ~Engine_Ext_CylinderMultiGrid(); + + void SetBarrier(boost::barrier* waitBase, boost::barrier* waitChild, boost::barrier* waitSync); + + virtual void DoPreVoltageUpdates(); + virtual void DoPostVoltageUpdates(); + virtual void Apply2Voltages(); + + virtual void DoPreCurrentUpdates(); + virtual void DoPostCurrentUpdates(); + virtual void Apply2Current(); + + virtual void SetEngine(Engine* eng); + +protected: + void SyncVoltages(); + void SyncCurrents(); + + Engine_CylinderMultiGrid* m_Eng_MG; + + boost::barrier *m_WaitOnBase; + boost::barrier *m_WaitOnChild; + boost::barrier *m_WaitOnSync; + + bool m_IsBase; +}; + +#endif // ENGINE_EXT_CYLINDERMULTIGRID_H diff --git a/FDTD/engine_multithread.h b/FDTD/engine_multithread.h index 07f9ec6..2b55eec 100644 --- a/FDTD/engine_multithread.h +++ b/FDTD/engine_multithread.h @@ -78,6 +78,7 @@ class Engine_Multithread : public Engine_SSE_Compressed { friend class NS_Engine_Multithread::thread; friend class NS_Engine_Multithread::thread_e_excitation; + friend class Engine_CylinderMultiGrid; public: static Engine_Multithread* New(const Operator_Multithread* op, unsigned int numThreads = 0); virtual ~Engine_Multithread(); diff --git a/FDTD/engine_sse.h b/FDTD/engine_sse.h index 7a53717..9389f84 100644 --- a/FDTD/engine_sse.h +++ b/FDTD/engine_sse.h @@ -52,6 +52,7 @@ protected: unsigned int numVectors; +public: //public access to the sse arrays for efficient extensions access... use careful... f4vector**** f4_volt; f4vector**** f4_curr; }; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 1e87285..4979183 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -467,7 +467,6 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo) { if (geo==NULL) return false; - Reset(); CSX = geo; CSRectGrid* grid=CSX->GetGrid(); @@ -584,7 +583,8 @@ int Operator::CalcECOperator() bool PEC[6]={1,1,1,1,1,1}; //exception for pml boundaries for (int n=0;n<6;++n) - PEC[n] = m_BC[n]!=3; + if ((m_BC[n]==3) || (m_BC[n]==-1)) + PEC[n] = false; ApplyElectricBC(PEC); InitExcitation(); @@ -620,8 +620,8 @@ void Operator::ApplyElectricBC(bool* dirs) GetVI(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n]; pos[n]=numLines[n]-1; - GetVV(n,pos[0],pos[1],pos[2]) = 0; // these are outside the FDTD-domain as defined by the main disc - GetVI(n,pos[0],pos[1],pos[2]) = 0; // these are outside the FDTD-domain as defined by the main disc + GetVV(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; // these are outside the FDTD-domain as defined by the main disc + GetVI(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; // these are outside the FDTD-domain as defined by the main disc GetVV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; GetVI(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; diff --git a/FDTD/operator.h b/FDTD/operator.h index b5cb728..56c2336 100644 --- a/FDTD/operator.h +++ b/FDTD/operator.h @@ -60,7 +60,7 @@ public: //! Set a forced timestep to use by the operator virtual void SetTimestep(double ts) {dT = ts;} double GetTimestep() const {return dT;}; - double GetNumberCells() const; + virtual double GetNumberCells() const; //! Returns the number of lines as needed for post-processing etc. (for the engine, use GetOriginalNumLines()) virtual unsigned int GetNumberOfLines(int ny) const {return numLines[ny];} diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index 8e7e588..5d4a59b 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -27,6 +27,7 @@ all special cases e.g. a closed alpha mesh or an included r=0 case is treated by */ class Operator_Cylinder : public Operator_Multithread { + friend class Operator_CylinderMultiGrid; public: static Operator_Cylinder* New(unsigned int numThreads = 0); virtual ~Operator_Cylinder(); diff --git a/FDTD/operator_cylindermultigrid.cpp b/FDTD/operator_cylindermultigrid.cpp new file mode 100644 index 0000000..f418ba4 --- /dev/null +++ b/FDTD/operator_cylindermultigrid.cpp @@ -0,0 +1,189 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY{} without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#include "operator_cylindermultigrid.h" +#include "engine_cylindermultigrid.h" +#include "operator_ext_cylinder.h" + +Operator_CylinderMultiGrid::Operator_CylinderMultiGrid(vector Split_Radii) : Operator_Cylinder() +{ + m_Split_Rad = Split_Radii.back(); + Split_Radii.pop_back(); +} + +Operator_CylinderMultiGrid::~Operator_CylinderMultiGrid() +{ +} + +Operator_CylinderMultiGrid* Operator_CylinderMultiGrid::New(vector Split_Radii, unsigned int numThreads) +{ + cout << "Create cylindrical multi grid FDTD operator " << endl; + Operator_CylinderMultiGrid* op = new Operator_CylinderMultiGrid(Split_Radii); + op->setNumThreads(numThreads); + op->Init(); + + return op; +} + +Engine* Operator_CylinderMultiGrid::CreateEngine() const +{ + Engine_CylinderMultiGrid* eng = Engine_CylinderMultiGrid::New(this,m_numThreads); + return eng; +} + +double Operator_CylinderMultiGrid::GetNumberCells() const +{ + if (numLines) + return (numLines[0]-m_Split_Pos)*(numLines[1])*(numLines[2]) + m_InnerOp->GetNumberCells(); + return 0; +} + +bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo) +{ + if (Operator_Cylinder::SetGeometryCSX(geo)==false) + return false; + + if (numLines[1]%2 != 1) + { + cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, number of line in alpha direction must be odd... found: " << numLines[1] << endl; + } + + m_Split_Pos = 0; + for (unsigned int n=0;nnumLines[0]-4)) + { + cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, split invalid..." << endl; + return false; + } + + CSRectGrid* grid = geo->GetGrid(); + + grid->ClearLines(0); + grid->ClearLines(1); + for (unsigned int n=0; nAddDiscLine(0,discLines[0][n]); + for (unsigned int n=0; nAddDiscLine(1,discLines[1][n]); + + if (m_InnerOp->SetGeometryCSX(CSX)==false) + return false; + + //restore grid to original mesh + grid->ClearLines(0); + grid->ClearLines(1); + for (unsigned int n=0; nAddDiscLine(0,discLines[0][n]); + for (unsigned int n=0; nAddDiscLine(1,discLines[1][n]); + + return true; +} + +void Operator_CylinderMultiGrid::Init() +{ + Operator_Cylinder::Init(); + m_InnerOp = Operator_Cylinder::New(m_numThreads); +} + +void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, vector &start, vector &stop) const +{ + if (numLines[0]= (numLines[0] - m_Split_Pos + 1)) + --numThreads; + + start.resize(numThreads); + stop.resize(numThreads); + + for (unsigned int n=0; nSetTimestep(dT); + + //calc inner child first + m_InnerOp->CalcECOperator(); + + dT = m_InnerOp->GetTimestep(); + + return Operator_Cylinder::CalcECOperator(); +} + +bool Operator_CylinderMultiGrid::SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) +{ + if (!m_InnerOp->Exc->setupExcitation(Excite,maxTS)) + return false; + return Exc->setupExcitation(Excite,maxTS); +} + +void Operator_CylinderMultiGrid::Reset() +{ + Operator_Cylinder::Reset(); + m_InnerOp->Reset(); +} + +void Operator_CylinderMultiGrid::SetBoundaryCondition(int* BCs) +{ + Operator_Cylinder::SetBoundaryCondition(BCs); + BCs[1] = 0; //always PEC in +r-direction + m_InnerOp->SetBoundaryCondition(BCs); +} + +void Operator_CylinderMultiGrid::AddExtension(Operator_Extension* op_ext) +{ + if (dynamic_cast(op_ext)) + { + Operator_Cylinder::AddExtension(op_ext); + return; + } + Operator_Extension* child_Ext = op_ext->Clone(m_InnerOp); + if (child_Ext==NULL) + { + cerr << "Operator_CylinderMultiGrid::AddExtension: Warning, extension: " << op_ext->GetExtensionName() << " can not be cloned for the child operator. Skipping Extension... " << endl; + return; + } + Operator_Cylinder::AddExtension(op_ext); + //give the copy to child + m_InnerOp->AddExtension(child_Ext); +} + +void Operator_CylinderMultiGrid::ShowStat() const +{ + m_InnerOp->ShowStat(); + m_InnerOp->ShowExtStat(); + Operator_Cylinder::ShowStat(); +} diff --git a/FDTD/operator_cylindermultigrid.h b/FDTD/operator_cylindermultigrid.h new file mode 100644 index 0000000..57612f3 --- /dev/null +++ b/FDTD/operator_cylindermultigrid.h @@ -0,0 +1,70 @@ +/* +* Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#ifndef OPERATOR_CYLINDERMULTIGRID_H +#define OPERATOR_CYLINDERMULTIGRID_H + +#include "operator_cylinder.h" + +//! This is a cylindrical FDTD operator using a simple multi-grid approach. +/*! + This cylindrical multi-grid operator itself is not calculating any real operator, instead it is hosting two separate "child" operators of type "Operator_Cylinder". + This operator class (or the corresponding engine) will perform the interpolation and connection between these two child-operator/engines. + One of the child operators itself may be another multi-grid operator to allow for a cascaded multi-grid approach. + */ +class Operator_CylinderMultiGrid : public Operator_Cylinder +{ + friend class Engine_CylinderMultiGrid; +public: + static Operator_CylinderMultiGrid* New(vector Split_Radii, unsigned int numThreads = 0); + virtual ~Operator_CylinderMultiGrid(); + + virtual double GetNumberCells() const; + + virtual Engine* CreateEngine() const; + + virtual bool SetGeometryCSX(ContinuousStructure* geo); + + virtual unsigned int GetSplitPos() const {return m_Split_Pos;} + + virtual int CalcECOperator(); + + virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS); + + virtual void SetBoundaryCondition(int* BCs); + + virtual void AddExtension(Operator_Extension* op_ext); + + Operator_Cylinder* GetInnerOperator() const {return m_InnerOp;} + + virtual void ShowStat() const; + +protected: + Operator_CylinderMultiGrid(vector Split_Radii); + virtual void Init(); +// virtual void InitOperator(); + virtual void Reset(); + + double m_Split_Rad; + unsigned int m_Split_Pos; + + Operator_Cylinder* m_InnerOp; + + virtual void CalcStartStopLines(unsigned int &numThreads, vector &start, vector &stop) const; +}; + +#endif // OPERATOR_CYLINDERMULTIGRID_H diff --git a/FDTD/operator_multithread.cpp b/FDTD/operator_multithread.cpp index 56f12eb..e0f1a16 100644 --- a/FDTD/operator_multithread.cpp +++ b/FDTD/operator_multithread.cpp @@ -76,6 +76,9 @@ void Operator_Multithread::Reset() void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector &start, vector &stop) const { + if (numLines[0]= numLines[0]) --numThreads; diff --git a/openEMS.pro b/openEMS.pro index a219e72..e26b3a1 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -81,7 +81,10 @@ SOURCES += main.cpp \ FDTD/engine_ext_lorentzmaterial.cpp \ FDTD/operator_ext_pml_sf.cpp \ FDTD/engine_ext_pml_sf.cpp \ - FDTD/processmodematch.cpp + FDTD/processmodematch.cpp \ + FDTD/operator_cylindermultigrid.cpp \ + FDTD/engine_cylindermultigrid.cpp \ + FDTD/engine_ext_cylindermultigrid.cpp HEADERS += tools/ErrorMsg.h \ tools/AdrOp.h \ tools/constants.h \ @@ -119,7 +122,10 @@ HEADERS += tools/ErrorMsg.h \ FDTD/engine_ext_lorentzmaterial.h \ FDTD/operator_ext_pml_sf.h \ FDTD/engine_ext_pml_sf.h \ - FDTD/processmodematch.h + FDTD/processmodematch.h \ + FDTD/operator_cylindermultigrid.h \ + FDTD/engine_cylindermultigrid.h \ + FDTD/engine_ext_cylindermultigrid.h QMAKE_CXXFLAGS_RELEASE = -O3 \ -g \ -march=native diff --git a/openems.cpp b/openems.cpp index bf69279..96707a6 100644 --- a/openems.cpp +++ b/openems.cpp @@ -19,6 +19,7 @@ #include #include "tools/array_ops.h" #include "FDTD/operator_cylinder.h" +#include "FDTD/operator_cylindermultigrid.h" #include "FDTD/engine_multithread.h" #include "FDTD/operator_multithread.h" #include "FDTD/operator_ext_mur_abc.h" @@ -293,7 +294,6 @@ int openEMS::SetupFDTD(const char* file) double maxTime=0; FDTD_Opts->QueryDoubleAttribute("MaxTime",&maxTime); - TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond"); if (BC==NULL) { @@ -319,7 +319,15 @@ int openEMS::SetupFDTD(const char* file) //*************** setup operator ************// if (CylinderCoords) { - FDTD_Op = Operator_Cylinder::New(m_engine_numThreads); + const char* radii = FDTD_Opts->Attribute("MultiGrid"); + if (radii) + { + string rad(radii); + FDTD_Op = Operator_CylinderMultiGrid::New(SplitString2Double(rad,','),m_engine_numThreads); + } + else + FDTD_Op = Operator_Cylinder::New(m_engine_numThreads); + CSX.SetCoordInputType(1); //tell CSX to use cylinder-coords } else if (m_engine == EngineType_SSE) { From 20ade0f05351091c883630819ffc05c32f0f8430 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Wed, 8 Sep 2010 16:07:28 +0200 Subject: [PATCH 13/13] new: enable cascaded multi-grids... incl. an example --- FDTD/engine.cpp | 1 - FDTD/engine_cylindermultigrid.cpp | 38 ++++--- FDTD/engine_cylindermultigrid.h | 6 +- FDTD/operator.cpp | 6 + FDTD/operator.h | 2 + FDTD/operator_cylindermultigrid.cpp | 38 ++++++- FDTD/operator_cylindermultigrid.h | 1 + matlab/examples/Coax_Cylindrical_MG.m | 155 ++++++++++++++++++++++++++ 8 files changed, 223 insertions(+), 24 deletions(-) create mode 100644 matlab/examples/Coax_Cylindrical_MG.m diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index 0d1e511..fc2a37e 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -48,7 +48,6 @@ Engine::~Engine() void Engine::Init() { - Reset(); numTS = 0; volt = Create_N_3DArray(numLines); curr = Create_N_3DArray(numLines); diff --git a/FDTD/engine_cylindermultigrid.cpp b/FDTD/engine_cylindermultigrid.cpp index c517f4e..ae94b76 100644 --- a/FDTD/engine_cylindermultigrid.cpp +++ b/FDTD/engine_cylindermultigrid.cpp @@ -36,15 +36,27 @@ Engine_CylinderMultiGrid::Engine_CylinderMultiGrid(const Operator_CylinderMultiG m_WaitOnChild = new boost::barrier(2); m_WaitOnSync = new boost::barrier(2); - Engine* eng = op->GetInnerOperator()->CreateEngine(); - m_InnerEngine = dynamic_cast(eng); - m_Eng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,true); m_Eng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync); m_Eng_Ext_MG->SetEngine(this); - m_InnerEng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,false); + + Engine* eng = op->GetInnerOperator()->CreateEngine(); + m_InnerEngine = dynamic_cast(eng); + + Engine_Ext_CylinderMultiGrid* m_InnerEng_Ext_MG = new Engine_Ext_CylinderMultiGrid(NULL,false); m_InnerEng_Ext_MG->SetBarrier(m_WaitOnBase, m_WaitOnChild, m_WaitOnSync); + // if already has a base extension, switch places ... seems to be faster... + for (size_t n=0;nm_Eng_exts.size();++n) + { + Engine_Ext_CylinderMultiGrid* eng_mg = dynamic_cast(m_InnerEngine->m_Eng_exts.at(n)); + if (eng_mg) + { + m_InnerEngine->m_Eng_exts.at(n) = m_InnerEng_Ext_MG; + m_InnerEng_Ext_MG = eng_mg; + break; + } + } m_InnerEngine->m_Eng_exts.push_back(m_InnerEng_Ext_MG); } @@ -83,10 +95,10 @@ void Engine_CylinderMultiGrid::Init() boost::thread *t = NULL; - t = new boost::thread( Engine_CylinderMultiGrid_Thread(this,m_startBarrier,m_stopBarrier,&m_Thread_NumTS) ); + t = new boost::thread( Engine_CylinderMultiGrid_Thread(this,m_startBarrier,m_stopBarrier,&m_Thread_NumTS, true) ); m_IteratorThread_Group.add_thread( t ); - t = new boost::thread( Engine_CylinderMultiGrid_Thread(m_InnerEngine,m_startBarrier,m_stopBarrier,&m_Thread_NumTS) ); + t = new boost::thread( Engine_CylinderMultiGrid_Thread(m_InnerEngine,m_startBarrier,m_stopBarrier,&m_Thread_NumTS, false) ); m_IteratorThread_Group.add_thread( t ); } @@ -107,12 +119,6 @@ bool Engine_CylinderMultiGrid::IterateTS(unsigned int iterTS) return true; } -void Engine_CylinderMultiGrid::InitExtensions() -{ - m_InnerEngine->InitExtensions(); - Engine_Multithread::InitExtensions(); -} - void Engine_CylinderMultiGrid::InterpolVoltChild2Base(unsigned int rzPlane) { //interpolate voltages from child engine to the base engine... @@ -229,11 +235,12 @@ void Engine_CylinderMultiGrid::InterpolCurrChild2Base(unsigned int rzPlane) } /****************************************************************************************/ -Engine_CylinderMultiGrid_Thread::Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS) +Engine_CylinderMultiGrid_Thread::Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS, bool isBase) { m_startBarrier = start; m_stopBarrier = stop; m_Eng=engine; + m_isBase=isBase; m_numTS = numTS; } @@ -243,7 +250,10 @@ void Engine_CylinderMultiGrid_Thread::operator()() while(*m_numTS>0) //m_numTS==0 request to terminate this thread... { - m_Eng->Engine_Multithread::IterateTS(*m_numTS); + if (m_isBase) + m_Eng->Engine_Multithread::IterateTS(*m_numTS); + else + m_Eng->IterateTS(*m_numTS); m_stopBarrier->wait(); //sync all workers after iterations are performed m_startBarrier->wait(); //wait for Base engine to start the iterations again ... } diff --git a/FDTD/engine_cylindermultigrid.h b/FDTD/engine_cylindermultigrid.h index 366f7b4..ac2f276 100644 --- a/FDTD/engine_cylindermultigrid.h +++ b/FDTD/engine_cylindermultigrid.h @@ -41,8 +41,6 @@ public: //! Iterate \a iterTS number of timesteps virtual bool IterateTS(unsigned int iterTS); - virtual void InitExtensions(); - protected: Engine_CylinderMultiGrid(const Operator_CylinderMultiGrid* op); const Operator_CylinderMultiGrid* Op_CMG; @@ -62,18 +60,18 @@ protected: boost::barrier *m_WaitOnSync; Engine_Ext_CylinderMultiGrid* m_Eng_Ext_MG; - Engine_Ext_CylinderMultiGrid* m_InnerEng_Ext_MG; }; class Engine_CylinderMultiGrid_Thread { public: - Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS); + Engine_CylinderMultiGrid_Thread( Engine_Multithread* engine, boost::barrier *start, boost::barrier *stop, volatile unsigned int* numTS, bool isBase); void operator()(); protected: Engine_Multithread *m_Eng; + bool m_isBase; boost::barrier *m_startBarrier; boost::barrier *m_stopBarrier; volatile unsigned int *m_numTS; diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 4979183..281ad8a 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -36,6 +36,7 @@ Operator::Operator() m_MeshType = ProcessFields::CARTESIAN_MESH; Exc = 0; dT = 0; + m_InvaildTimestep = false; } Operator::~Operator() @@ -535,6 +536,7 @@ int Operator::CalcECOperator() if (Calc_EC()==0) return -1; + m_InvaildTimestep = false; opt_dT = 0; if (dT>0) { @@ -542,7 +544,11 @@ int Operator::CalcECOperator() CalcTimestep(); opt_dT = dT; if (dT Split_Radii) : Operator_Cylinder() { - m_Split_Rad = Split_Radii.back(); - Split_Radii.pop_back(); + m_Split_Radii = Split_Radii; + m_Split_Rad = m_Split_Radii.back(); + m_Split_Radii.pop_back(); } Operator_CylinderMultiGrid::~Operator_CylinderMultiGrid() @@ -60,6 +61,18 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo) if (numLines[1]%2 != 1) { cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, number of line in alpha direction must be odd... found: " << numLines[1] << endl; + exit(0); + } + + //check if mesh is homogenous in alpha-direction + double diff=discLines[1][1]-discLines[1][0]; + for (unsigned int n=2;n 1e-10) + { + cerr << "Operator_CylinderMultiGrid::SetGeometryCSX: Error, mesh has to be homogenous in alpha direction for multi grid engine, violation found at: " << n << endl; + exit(0); + } } m_Split_Pos = 0; @@ -105,7 +118,11 @@ bool Operator_CylinderMultiGrid::SetGeometryCSX(ContinuousStructure* geo) void Operator_CylinderMultiGrid::Init() { Operator_Cylinder::Init(); - m_InnerOp = Operator_Cylinder::New(m_numThreads); + + if (m_Split_Radii.empty()) + m_InnerOp = Operator_Cylinder::New(m_numThreads); + else + m_InnerOp = Operator_CylinderMultiGrid::New(m_Split_Radii,m_numThreads); } void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, vector &start, vector &stop) const @@ -132,6 +149,7 @@ void Operator_CylinderMultiGrid::CalcStartStopLines(unsigned int &numThreads, ve int Operator_CylinderMultiGrid::CalcECOperator() { + int retCode=0; if (dT) m_InnerOp->SetTimestep(dT); @@ -140,12 +158,22 @@ int Operator_CylinderMultiGrid::CalcECOperator() dT = m_InnerOp->GetTimestep(); - return Operator_Cylinder::CalcECOperator(); + retCode = Operator_Cylinder::CalcECOperator(); + if (GetTimestepValid()==false) + { + cerr << "Operator_CylinderMultiGrid::CalcECOperator(): Warning, timestep invalid... resetting..." << endl; + dT = opt_dT; + m_InnerOp->SetTimestep(dT); + m_InnerOp->CalcECOperator(); + return Operator_Cylinder::CalcECOperator(); + } + + return retCode; } bool Operator_CylinderMultiGrid::SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) { - if (!m_InnerOp->Exc->setupExcitation(Excite,maxTS)) + if (!m_InnerOp->SetupExcitation(Excite,maxTS)) return false; return Exc->setupExcitation(Excite,maxTS); } diff --git a/FDTD/operator_cylindermultigrid.h b/FDTD/operator_cylindermultigrid.h index 57612f3..6fb7503 100644 --- a/FDTD/operator_cylindermultigrid.h +++ b/FDTD/operator_cylindermultigrid.h @@ -60,6 +60,7 @@ protected: virtual void Reset(); double m_Split_Rad; + vector m_Split_Radii; unsigned int m_Split_Pos; Operator_Cylinder* m_InnerOp; diff --git a/matlab/examples/Coax_Cylindrical_MG.m b/matlab/examples/Coax_Cylindrical_MG.m new file mode 100644 index 0000000..84a1668 --- /dev/null +++ b/matlab/examples/Coax_Cylindrical_MG.m @@ -0,0 +1,155 @@ +close all +clear +clc + +%example for an cylindrical mesh, modeling a coaxial cable +% this example is using a multi-grid approach + + +%% setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Settings = []; +Settings.LogFile = 'openEMS.log'; + +physical_constants + +f0 = 0.5e9; +epsR = 1; %material filling + +length = 1000; +port_dist = length/2; +rad_i = 10; %inner radius +rad_a = 200; %outer radius +partial = 0.5; %e.g. 0.5 means only one half of a coax, should be <1 or change boundary cond. +max_mesh = 10 / sqrt(epsR); +max_alpha = max_mesh; +N_alpha = ceil(rad_a * 2*pi * partial / max_alpha); + +%make it even... +N_alpha = N_alpha + mod(N_alpha,2); +%make sure it is multiple of 4, needed for 2 multi-grid steps +N_alpha = ceil((N_alpha)/4) *4 + 1; + +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --numThreads=1']; + +def_refSimu = 0; % do a reference simulation without the multi-grid + +%% setup done %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +if (def_refSimu>0) + Sim_Path = 'tmp_ref'; +else + Sim_Path = 'tmp'; +end +Sim_CSX = 'coax.xml'; + +if (exist(Sim_Path,'dir')) + rmdir(Sim_Path,'s'); +end +mkdir(Sim_Path); + +%setup FDTD parameter +if (def_refSimu>0) + FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',5 ); +else + FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',5 ,'MultiGrid','60,120'); +end +FDTD = SetGaussExcite(FDTD,f0,f0); +BC = [0 0 1 1 2 2]; +FDTD = SetBoundaryCond(FDTD,BC); + +mesh_res = [max_mesh 2*pi*partial/N_alpha max_mesh]; + +%setup CSXCAD geometry +CSX = InitCSX(); +mesh.x = SmoothMeshLines([rad_i rad_a],mesh_res(1)); +mesh.y = linspace(-pi*partial,pi*partial,N_alpha); +mesh.z = SmoothMeshLines([0 port_dist length],mesh_res(3)); +CSX = DefineRectGrid(CSX, 1e-3,mesh); + +start = [rad_i mesh.y(1) mesh.z(3)]; +stop = [rad_a mesh.y(end) mesh.z(3)]; + +CSX = AddExcitation(CSX,'excite',0,[1 0 0]); +weight{1} = '1/rho'; +weight{2} = 0; +weight{3} = 0; +CSX = SetExcitationWeight(CSX, 'excite', weight ); +CSX = AddBox(CSX,'excite',0 ,start,stop); + + +start = [mesh.x(1) mesh.y(1) mesh.z(1)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end)]; +CSX = AddMaterial(CSX,'material'); +CSX = SetMaterialProperty(CSX,'material','Epsilon',epsR); +CSX = AddBox(CSX,'material',0 ,start,stop); + +%dump +CSX = AddDump(CSX,'Et_rz_','DumpMode',0); +start = [mesh.x(1) 0 mesh.z(1)]; +stop = [mesh.x(end) 0 mesh.z(end)]; +CSX = AddBox(CSX,'Et_rz_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_rz_','DumpType',1,'DumpMode',0); +CSX = AddBox(CSX,'Ht_rz_',0 , start,stop); + +CSX = AddDump(CSX,'Et_','DumpType',0,'DumpMode',0); +start = [mesh.x(1) mesh.y(1) length/2]; +stop = [mesh.x(end) mesh.y(end) length/2]; +CSX = AddBox(CSX,'Et_',0,start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); +start = [mesh.x(1) mesh.y(1) length/2]; +stop = [mesh.x(end) mesh.y(end) length/2]; +CSX = AddBox(CSX,'Ht_',0,start,stop); + +% voltage calc (take a voltage average to be at the same spot as the +% current calculation) +CSX = AddProbe(CSX,'ut1_1',0); +start = [ rad_i 0 port_dist ];stop = [ rad_a 0 port_dist ]; +CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); +CSX = AddProbe(CSX,'ut1_2',0); +start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ]; +CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); + +% current calc +CSX = AddProbe(CSX,'it1',1); +mid = 75; +start = [ 0 mesh.y(1) port_dist+mesh_res(3)/2 ];stop = [ mid mesh.y(end) port_dist+mesh_res(3)/2 ]; +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts, Settings) + +%% +close all +freq = linspace(0,2*f0,201); +UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path,freq); +u_f = (UI.FD{1}.val + UI.FD{2}.val)/2; %averaging voltages to fit current +i_f = UI.FD{3}.val / partial; + +% plot(UI.TD{1}.t,UI.TD{1}.val); +% grid on; +% +% figure +% plot(UI.TD{3}.t,UI.TD{3}.val); +% grid on; + +%plot Z_L compare +figure +ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax +plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g','Linewidth',3); +hold on; +grid on; +Z = u_f./i_f; +plot(UI.FD{1}.f,real(Z),'k--','Linewidth',2); +plot(UI.FD{1}.f,imag(Z),'r-','Linewidth',2); +xlim([0 2*f0]); +legend('Z_L - analytic','\Re\{Z\} - FDTD','\Im\{Z\} - FDTD','Location','Best'); + +