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) {