diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp
index 5b25ba2..09c011a 100644
--- a/FDTD/engine.cpp
+++ b/FDTD/engine.cpp
@@ -18,10 +18,18 @@
#include "engine.h"
#include "tools/array_ops.h"
-Engine::Engine(Operator* op)
+//! \brief construct an Engine instance
+//! it's the responsibility of the caller to free the returned pointer
+Engine* Engine::createEngine(const Operator* op)
+{
+ Engine* e = new Engine(op);
+ e->Init();
+ return e;
+}
+
+Engine::Engine(const Operator* op)
{
Op = op;
- Init();
}
Engine::~Engine()
diff --git a/FDTD/engine.h b/FDTD/engine.h
index e1660b7..7068ac2 100644
--- a/FDTD/engine.h
+++ b/FDTD/engine.h
@@ -23,7 +23,7 @@
class Engine
{
public:
- Engine(Operator* op);
+ static Engine* createEngine(const Operator* op);
virtual ~Engine();
virtual void Init();
@@ -32,13 +32,14 @@ public:
//!Iterate a number of timesteps
virtual bool IterateTS(unsigned int iterTS);
- unsigned int GetNumberOfTimesteps() {return numTS;};
+ virtual unsigned int GetNumberOfTimesteps() {return numTS;};
virtual FDTD_FLOAT**** GetVoltages() {return volt;};
virtual FDTD_FLOAT**** GetCurrents() {return curr;};
protected:
- Operator* Op;
+ Engine(const Operator* op);
+ const Operator* Op;
FDTD_FLOAT**** volt;
FDTD_FLOAT**** curr;
diff --git a/FDTD/engine_multithread.cpp b/FDTD/engine_multithread.cpp
new file mode 100644
index 0000000..42ab92a
--- /dev/null
+++ b/FDTD/engine_multithread.cpp
@@ -0,0 +1,301 @@
+/*
+* Copyright (C) 2010 Sebastian Held (sebastian.held@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 .
+*/
+
+//#define ENABLE_DEBUG_TIME
+
+#ifdef ENABLE_DEBUG_TIME
+ #define DEBUG_TIME(x) x;
+#else
+ #define DEBUG_TIME(x) ;
+#endif
+
+
+
+#include "engine_multithread.h"
+#include "tools/array_ops.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+#include "boost/date_time/gregorian/gregorian.hpp"
+#include
+
+//! \brief construct an Engine_Multithread instance
+//! it's the responsibility of the caller to free the returned pointer
+Engine_Multithread* Engine_Multithread::createEngine(const Operator* op, unsigned int numThreads)
+{
+ Engine_Multithread* e = new Engine_Multithread(op);
+ e->setNumThreads( numThreads );
+ e->Init();
+ return e;
+}
+
+Engine_Multithread::Engine_Multithread(const Operator* op) : Engine(op)
+{
+}
+
+Engine_Multithread::~Engine_Multithread()
+{
+#ifdef ENABLE_DEBUG_TIME
+ NS_Engine_Multithread::DBG().cout() << "Engine_Multithread::~Engine_Multithread()" << endl;
+ std::map >::iterator it;
+ for (it=m_timer_list.begin(); it!=m_timer_list.end(); it++) {
+ NS_Engine_Multithread::DBG().cout() << "*** DEBUG Thread: " << it->first << std::endl;
+ std::vector::iterator it2;
+ for (it2=it->second.begin(); it2second.end();) {
+ NS_Engine_Multithread::DBG().cout() << "after voltage update, before barrier1: " << fixed << setprecision(6) << *(it2++) << std::endl;
+ NS_Engine_Multithread::DBG().cout() << "after barrier1, before barrier2: " << fixed << setprecision(6) << *(it2++) << std::endl;
+ NS_Engine_Multithread::DBG().cout() << "after barrier2, before current update: " << fixed << setprecision(6) << *(it2++) << std::endl;
+ NS_Engine_Multithread::DBG().cout() << "after current update, before barrier3: " << fixed << setprecision(6) << *(it2++) << std::endl;
+ NS_Engine_Multithread::DBG().cout() << "after barrier3: " << fixed << setprecision(6) << *(it2++) << std::endl;
+ }
+ }
+#endif
+
+ Reset();
+}
+
+void Engine_Multithread::setNumThreads( unsigned int numThreads )
+{
+ m_numThreads = numThreads;
+}
+
+void Engine_Multithread::Init()
+{
+ Engine::Init(); // gets cleaned up by Engine::~Engine()
+
+ // initialize threads
+ m_stopThreads = false;
+ if (m_numThreads == 0)
+ m_numThreads = boost::thread::hardware_concurrency();
+ cout << "using " << m_numThreads << " threads" << std::endl;
+ m_barrier1 = new boost::barrier(m_numThreads+1); // numThread workers + 1 excitation thread
+ m_barrier2 = new boost::barrier(m_numThreads+1); // numThread workers + 1 excitation thread
+ m_barrier3 = new boost::barrier(m_numThreads); // numThread workers
+ 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)Op->numLines[0] / (float)m_numThreads);
+ for (unsigned int n=0; nnumLines[0]-1;
+ stop_h = stop-1;
+ }
+ //NS_Engine_Multithread::DBG().cout() << "###DEBUG## Thread " << n << ": start=" << start << " stop=" << stop << " stop_h=" << stop_h << std::endl;
+ boost::thread *t = new boost::thread( NS_Engine_Multithread::thread(this,start,stop,stop_h,n) );
+ m_thread_group.add_thread( t );
+ }
+ boost::thread *t = new boost::thread( NS_Engine_Multithread::thread_e_excitation(this) );
+ m_thread_group.add_thread( t );
+}
+
+void Engine_Multithread::Reset()
+{
+ if (!m_stopThreads) {
+ // prevent multiple invocations
+
+ // stop the threads
+ //NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl;
+ m_iterTS = 1;
+ m_startBarrier->wait(); // start the threads
+ m_stopThreads = true;
+ m_stopBarrier->wait(); // wait for the threads to finish
+ m_thread_group.join_all(); // wait for termination
+ delete m_barrier1; m_barrier1 = 0;
+ delete m_barrier2; m_barrier2 = 0;
+ delete m_barrier3; m_barrier3 = 0;
+ delete m_startBarrier; m_startBarrier = 0;
+ delete m_stopBarrier; m_stopBarrier = 0;
+ }
+
+ Engine::Reset();
+}
+
+bool Engine_Multithread::IterateTS(unsigned int iterTS)
+{
+ m_iterTS = iterTS;
+
+ //cout << "bool Engine_Multithread::IterateTS(): starting threads ...";
+ m_startBarrier->wait(); // start the threads
+
+ //cout << "... threads started";
+
+ m_stopBarrier->wait(); // wait for the threads to finish time steps
+ return true;
+}
+
+//
+// *************************************************************************************************************************
+//
+namespace NS_Engine_Multithread {
+
+thread::thread( Engine_Multithread* ptr, unsigned int start, unsigned int stop, unsigned int stop_h, unsigned int threadID )
+{
+ m_enginePtr = ptr;
+ m_start = start;
+ m_stop = stop;
+ m_stop_h = stop_h;
+ m_threadID = threadID;
+}
+
+void thread::operator()()
+{
+ //std::cout << "thread::operator() Parameters: " << m_start << " " << m_stop << std::endl;
+ //DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") started." << endl;
+
+ unsigned int pos[3];
+ bool shift[3];
+
+ while (!m_enginePtr->m_stopThreads) {
+ // wait for start
+ //DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") waiting..." << endl;
+ m_enginePtr->m_startBarrier->wait();
+ //cout << "Thread " << boost::this_thread::get_id() << " waiting... started." << endl;
+
+ DEBUG_TIME( Timer timer1 );
+
+ for (unsigned int iter=0;iterm_iterTS;++iter)
+ {
+ //voltage updates
+ for (pos[0]=m_start;pos[0]<=m_stop;++pos[0])
+ {
+ shift[0]=pos[0];
+ for (pos[1]=0;pos[1]Op->numLines[1];++pos[1])
+ {
+ shift[1]=pos[1];
+ for (pos[2]=0;pos[2]Op->numLines[2];++pos[2])
+ {
+ shift[2]=pos[2];
+ //do the updates here
+ //for x
+ m_enginePtr->volt[0][pos[0]][pos[1]][pos[2]] *= m_enginePtr->Op->vv[0][pos[0]][pos[1]][pos[2]];
+ m_enginePtr->volt[0][pos[0]][pos[1]][pos[2]] += m_enginePtr->Op->vi[0][pos[0]][pos[1]][pos[2]] * ( m_enginePtr->curr[2][pos[0]][pos[1]][pos[2]] - m_enginePtr->curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - m_enginePtr->curr[1][pos[0]][pos[1]][pos[2]] + m_enginePtr->curr[1][pos[0]][pos[1]][pos[2]-shift[2]]);
+
+ //for y
+ m_enginePtr->volt[1][pos[0]][pos[1]][pos[2]] *= m_enginePtr->Op->vv[1][pos[0]][pos[1]][pos[2]];
+ m_enginePtr->volt[1][pos[0]][pos[1]][pos[2]] += m_enginePtr->Op->vi[1][pos[0]][pos[1]][pos[2]] * ( m_enginePtr->curr[0][pos[0]][pos[1]][pos[2]] - m_enginePtr->curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - m_enginePtr->curr[2][pos[0]][pos[1]][pos[2]] + m_enginePtr->curr[2][pos[0]-shift[0]][pos[1]][pos[2]]);
+
+ //for x
+ m_enginePtr->volt[2][pos[0]][pos[1]][pos[2]] *= m_enginePtr->Op->vv[2][pos[0]][pos[1]][pos[2]];
+ m_enginePtr->volt[2][pos[0]][pos[1]][pos[2]] += m_enginePtr->Op->vi[2][pos[0]][pos[1]][pos[2]] * ( m_enginePtr->curr[1][pos[0]][pos[1]][pos[2]] - m_enginePtr->curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - m_enginePtr->curr[0][pos[0]][pos[1]][pos[2]] + m_enginePtr->curr[0][pos[0]][pos[1]-shift[1]][pos[2]]);
+ }
+ }
+ }
+
+ // record time
+ DEBUG_TIME( m_enginePtr->m_timer_list[boost::this_thread::get_id()].push_back( timer1.elapsed() ); )
+
+ //cout << "Thread " << boost::this_thread::get_id() << " m_barrier1 waiting..." << endl;
+ m_enginePtr->m_barrier1->wait();
+
+ // record time
+ DEBUG_TIME( m_enginePtr->m_timer_list[boost::this_thread::get_id()].push_back( timer1.elapsed() ); )
+
+ // e-field excitation (thread thread_e_excitation)
+
+ m_enginePtr->m_barrier2->wait();
+ // e_excitation finished
+
+ // record time
+ DEBUG_TIME( m_enginePtr->m_timer_list[boost::this_thread::get_id()].push_back( timer1.elapsed() ); )
+
+ //current updates
+ for (pos[0]=m_start;pos[0]<=m_stop_h;++pos[0])
+ {
+ for (pos[1]=0;pos[1]Op->numLines[1]-1;++pos[1])
+ {
+ for (pos[2]=0;pos[2]Op->numLines[2]-1;++pos[2])
+ {
+ //do the updates here
+ //for x
+ m_enginePtr->curr[0][pos[0]][pos[1]][pos[2]] *= m_enginePtr->Op->ii[0][pos[0]][pos[1]][pos[2]];
+ m_enginePtr->curr[0][pos[0]][pos[1]][pos[2]] += m_enginePtr->Op->iv[0][pos[0]][pos[1]][pos[2]] * ( m_enginePtr->volt[2][pos[0]][pos[1]][pos[2]] - m_enginePtr->volt[2][pos[0]][pos[1]+1][pos[2]] - m_enginePtr->volt[1][pos[0]][pos[1]][pos[2]] + m_enginePtr->volt[1][pos[0]][pos[1]][pos[2]+1]);
+
+ //for y
+ m_enginePtr->curr[1][pos[0]][pos[1]][pos[2]] *= m_enginePtr->Op->ii[1][pos[0]][pos[1]][pos[2]];
+ m_enginePtr->curr[1][pos[0]][pos[1]][pos[2]] += m_enginePtr->Op->iv[1][pos[0]][pos[1]][pos[2]] * ( m_enginePtr->volt[0][pos[0]][pos[1]][pos[2]] - m_enginePtr->volt[0][pos[0]][pos[1]][pos[2]+1] - m_enginePtr->volt[2][pos[0]][pos[1]][pos[2]] + m_enginePtr->volt[2][pos[0]+1][pos[1]][pos[2]]);
+
+ //for x
+ m_enginePtr->curr[2][pos[0]][pos[1]][pos[2]] *= m_enginePtr->Op->ii[2][pos[0]][pos[1]][pos[2]];
+ m_enginePtr->curr[2][pos[0]][pos[1]][pos[2]] += m_enginePtr->Op->iv[2][pos[0]][pos[1]][pos[2]] * ( m_enginePtr->volt[1][pos[0]][pos[1]][pos[2]] - m_enginePtr->volt[1][pos[0]+1][pos[1]][pos[2]] - m_enginePtr->volt[0][pos[0]][pos[1]][pos[2]] + m_enginePtr->volt[0][pos[0]][pos[1]+1][pos[2]]);
+ }
+ }
+ }
+
+ // record time
+ DEBUG_TIME( m_enginePtr->m_timer_list[boost::this_thread::get_id()].push_back( timer1.elapsed() ); )
+
+ m_enginePtr->m_barrier3->wait();
+
+ // record time
+ DEBUG_TIME( m_enginePtr->m_timer_list[boost::this_thread::get_id()].push_back( timer1.elapsed() ); )
+
+ //soft current excitation here (H-field excite)
+
+ if (m_threadID == 0)
+ ++m_enginePtr->numTS; // only the first thread increments numTS
+ }
+
+ m_enginePtr->m_stopBarrier->wait();
+ }
+
+ //DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") finished." << endl;
+}
+
+} // namespace
+
+//
+// *************************************************************************************************************************
+//
+namespace NS_Engine_Multithread {
+
+thread_e_excitation::thread_e_excitation( Engine_Multithread* ptr )
+{
+ m_enginePtr = ptr;
+}
+
+void thread_e_excitation::operator()()
+{
+ //std::cout << "thread_e_excitation::operator()" << std::endl;
+ //DBG().cout() << "Thread e_excitation (" << boost::this_thread::get_id() << ") started." << endl;
+
+ int exc_pos;
+ const unsigned int E_Exc_Count = m_enginePtr->Op->E_Exc_Count;
+
+ while (!m_enginePtr->m_stopThreads)
+ {
+ // waiting on NS_Engine_Multithread::thread
+ m_enginePtr->m_barrier1->wait();
+
+ // soft voltage excitation here (E-field excite)
+ for (unsigned int n=0;nnumTS - (int)m_enginePtr->Op->E_Exc_delay[n];
+ exc_pos*= (exc_pos>0 && exc_pos<=(int)m_enginePtr->Op->ExciteLength);
+ m_enginePtr->volt[m_enginePtr->Op->E_Exc_dir[n]][m_enginePtr->Op->E_Exc_index[0][n]][m_enginePtr->Op->E_Exc_index[1][n]][m_enginePtr->Op->E_Exc_index[2][n]] += m_enginePtr->Op->E_Exc_amp[n]*m_enginePtr->Op->ExciteSignal[exc_pos];
+ }
+
+ // continue NS_Engine_Multithread::thread
+ m_enginePtr->m_barrier2->wait();
+ }
+
+ //DBG().cout() << "Thread e_excitation (" << boost::this_thread::get_id() << ") finished." << endl;
+}
+
+} // namespace
diff --git a/FDTD/engine_multithread.h b/FDTD/engine_multithread.h
new file mode 100644
index 0000000..ef3623d
--- /dev/null
+++ b/FDTD/engine_multithread.h
@@ -0,0 +1,99 @@
+/*
+* Copyright (C) 2010 Sebastian Held (sebastian.held@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_MULTITHREAD_H
+#define ENGINE_MULTITHREAD_H
+
+#include "operator.h"
+#include "engine.h"
+
+#include
+#include
+#include
+#include
+
+class Engine_Multithread;
+
+namespace NS_Engine_Multithread {
+
+ class DBG { // debug
+ public:
+ DBG() {}
+ ~DBG() { std::cout << os.str();}
+ std::ostringstream& cout() {return os;}
+ protected:
+ std::ostringstream os;
+ };
+
+ class Timer { //debug
+ public:
+ Timer() {gettimeofday(&t1,NULL);}
+ double elapsed() {gettimeofday(&t2,NULL); return (t2.tv_sec-t1.tv_sec) + (t2.tv_usec-t1.tv_usec)*1e-6;}
+ protected:
+ timeval t1,t2;
+ };
+
+ class thread {
+ public:
+ thread( Engine_Multithread* ptr, unsigned int start, unsigned int stop, unsigned int stop_h, unsigned int threadID );
+ void operator()();
+
+ protected:
+ unsigned int m_start, m_stop, m_stop_h, m_threadID;
+ Engine_Multithread *m_enginePtr;
+ };
+
+ class thread_e_excitation {
+ public:
+ thread_e_excitation( Engine_Multithread* ptr);
+ void operator()();
+
+ protected:
+ Engine_Multithread *m_enginePtr;
+ };
+} // namespace
+
+
+class Engine_Multithread : public Engine
+{
+ friend class NS_Engine_Multithread::thread;
+ friend class NS_Engine_Multithread::thread_e_excitation;
+public:
+ static Engine_Multithread* createEngine(const Operator* op, unsigned int numThreads = 0);
+ virtual ~Engine_Multithread();
+
+ virtual void setNumThreads( unsigned int numThreads );
+ virtual void Init();
+ virtual void Reset();
+
+ //!Iterate a number of timesteps
+ virtual bool IterateTS(unsigned int iterTS);
+
+protected:
+ Engine_Multithread(const Operator* op);
+ boost::thread_group m_thread_group;
+ boost::barrier *m_barrier1, *m_barrier2, *m_barrier3, *m_startBarrier, *m_stopBarrier;
+ volatile unsigned int m_iterTS;
+ unsigned int m_numThreads; //!< number of worker threads
+ volatile bool m_stopThreads;
+
+#ifdef ENABLE_DEBUG_TIME
+ std::map > m_timer_list;
+#endif
+};
+
+#endif // ENGINE_MULTITHREAD_H
diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index a3fb64a..2fd1638 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -205,14 +205,14 @@ struct Operator::Grid_Path Operator::FindPath(double start[], double stop[])
return path;
}
-double Operator::GetNumberCells()
+double Operator::GetNumberCells() const
{
if (numLines)
return (numLines[0])*(numLines[1])*(numLines[2]); //it's more like number of nodes???
return 0;
}
-void Operator::ShowStat()
+void Operator::ShowStat() const
{
unsigned int OpSize = 12*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
unsigned int FieldSize = 6*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
@@ -248,6 +248,39 @@ unsigned int Operator::CalcGaussianPulsExcitation(double f0, double fc)
return GetNyquistNum(f0+fc);
}
+unsigned int Operator::CalcDiracPulsExcitation()
+{
+ if (dT==0) return 0;
+
+ ExciteLength = 1;
+// cerr << "Operator::CalcDiracPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
+ delete[] ExciteSignal;
+ ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
+ ExciteSignal[0]=0.0;
+ ExciteSignal[1]=1.0;
+
+ // FIXME GetNyquistNum() has side-effects!
+ m_nyquistTS = 1;
+
+ return m_nyquistTS;
+}
+
+unsigned int Operator::CalcStepExcitation()
+{
+ if (dT==0) return 0;
+
+ ExciteLength = 1;
+ delete[] ExciteSignal;
+ ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
+ ExciteSignal[0]=1.0;
+ ExciteSignal[1]=1.0;
+
+ // FIXME GetNyquistNum() has side-effects!
+ m_nyquistTS = 1;
+
+ return m_nyquistTS;
+}
+
unsigned int Operator::CalcSinusExcitation(double f0, int nTS)
{
if (dT==0) return 0;
diff --git a/FDTD/operator.h b/FDTD/operator.h
index 12ed41a..fab5c27 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -41,15 +41,19 @@ public:
virtual unsigned int CalcGaussianPulsExcitation(double f0, double fc);
//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps \return number of Nyquist timesteps
virtual unsigned int CalcSinusExcitation(double f0, int nTS);
+ //! Calculate a dirac impuls excitation \return number of Nyquist timesteps
+ virtual unsigned int CalcDiracPulsExcitation();
+ //! Calculate a step excitation \return number of Nyquist timesteps
+ virtual unsigned int CalcStepExcitation();
virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
virtual void ApplyMagneticBC(bool* dirs);
- double GetTimestep() {return dT;};
+ double GetTimestep() const {return dT;};
unsigned int GetNyquistNum(double fmax);
- double GetNumberCells();
+ double GetNumberCells() const;
- void ShowStat();
+ void ShowStat() const;
void DumpOperator2File(string filename);
void DumpMaterial2File(string filename);
diff --git a/FDTD/processing.h b/FDTD/processing.h
index a5cd4e3..f9dd017 100644
--- a/FDTD/processing.h
+++ b/FDTD/processing.h
@@ -31,7 +31,7 @@ public:
virtual void DefineStartStopCoord(double* dstart, double* dstop);
- void SetProcessInterval(unsigned int interval) {ProcessInterval=interval;}
+ void SetProcessInterval(unsigned int interval) {ProcessInterval=max((unsigned int)1,interval);}
virtual int Process() {return GetNextInterval();}
//! If Disabled Process() will do nothing...
diff --git a/TESTSUITE/combinedtests/Coax.m b/TESTSUITE/combinedtests/Coax.m
new file mode 100644
index 0000000..1ea2009
--- /dev/null
+++ b/TESTSUITE/combinedtests/Coax.m
@@ -0,0 +1,152 @@
+function pass = Coax
+
+physical_constants;
+
+
+ENABLE_PLOTS = 1;
+CLEANUP = 0; % if enabled and result is PASS, remove simulation folder
+STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
+
+% LIMITS
+upper_error = 0.036; % max +3.6%
+lower_error = 0; % max -0.0%
+
+% structure
+abs_length = 250;
+length = 1000;
+coax_rad_i = 100;
+coax_rad_ai = 230;
+coax_rad_aa = 240;
+mesh_res = [5 5 5];
+f_start = 0;
+f_stop = 1e9;
+
+Sim_Path = 'tmp';
+Sim_CSX = 'coax.xml';
+
+[status,message,messageid]=mkdir(Sim_Path);
+
+%setup FDTD parameter
+FDTD = InitFDTD(5e5,1e-6);
+FDTD = SetGaussExcite(FDTD,(f_stop-f_start)/2,(f_stop-f_start)/2);
+BC = [1 1 1 1 1 1] * 0;
+FDTD = SetBoundaryCond(FDTD,BC);
+
+%setup CSXCAD geometry
+CSX = InitCSX();
+mesh.x = -2.5*mesh_res(1)-coax_rad_aa : mesh_res(1) : coax_rad_aa+2.5*mesh_res(1);
+mesh.y = mesh.x;
+mesh.z = 0 : mesh_res(3) : length;
+CSX = DefineRectGrid(CSX, 1e-3,mesh);
+
+%create copper
+CSX = AddMetal(CSX,'PEC');
+
+%%%fake pml
+finalKappa = 0.3/abs_length^4;
+finalSigma = finalKappa*MUE0/EPS0;
+CSX = AddMaterial(CSX,'pml');
+CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa);
+CSX = SetMaterialProperty(CSX,'pml','Sigma',finalSigma);
+CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
+CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
+
+%%% coax
+start = [0, 0 , 0];stop = [0, 0 , length];
+CSX = AddCylinder(CSX,'PEC',0 ,start,stop,coax_rad_i); % inner conductor
+CSX = AddCylindricalShell(CSX,'PEC',0 ,start,stop,0.5*(coax_rad_aa+coax_rad_ai),(coax_rad_aa-coax_rad_ai)); % outer conductor
+
+%%% add PML
+start(3) = length-abs_length;
+CSX = AddCylindricalShell(CSX,'pml',0 ,start,stop,0.5*(coax_rad_i+coax_rad_ai),(coax_rad_ai-coax_rad_i));
+start(3) = 0; stop(3)=mesh_res(1)/2;
+CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
+weight{1} = '(x)/(x*x+y*y)';
+weight{2} = 'y/pow(rho,2)';
+weight{3} = 0;
+CSX = SetExcitationWeight(CSX, 'excite', weight );
+CSX = AddCylindricalShell(CSX,'excite',0 ,start,stop,0.5*(coax_rad_i+coax_rad_ai),(coax_rad_ai-coax_rad_i));
+
+%dump
+CSX = AddDump(CSX,'Et_',0,2);
+start = [mesh.x(1) , 0 , mesh.z(1)];
+stop = [mesh.x(end) , 0 , mesh.z(end)];
+CSX = AddBox(CSX,'Et_',0 , start,stop);
+
+CSX = AddDump(CSX,'Ht_',1,2);
+CSX = AddBox(CSX,'Ht_',0,start,stop);
+
+%voltage calc
+CSX = AddProbe(CSX,'ut1',0);
+start = [ coax_rad_i 0 length/2 ];stop = [ coax_rad_ai 0 length/2 ];
+CSX = AddBox(CSX,'ut1', 0 ,start,stop);
+
+%current calc
+CSX = AddProbe(CSX,'it1',1);
+mid = 0.5*(coax_rad_i+coax_rad_ai);
+start = [ -mid -mid length/2 ];stop = [ mid mid length/2 ];
+CSX = AddBox(CSX,'it1', 0 ,start,stop);
+
+%Write openEMS compatible xml-file
+WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
+
+%cd to working dir and run openEMS
+savePath = pwd();
+cd(Sim_Path); %cd to working dir
+invoke_openEMS( Sim_CSX );
+UI = ReadUI( {'ut1','it1'} );
+cd(savePath);
+
+
+
+%
+% analysis
+%
+
+f = UI.FD{2}.f;
+u = UI.FD{1}.val;
+i = UI.FD{2}.val;
+
+f_idx_start = interp1( f, 1:numel(f), f_start, 'nearest' );
+f_idx_stop = interp1( f, 1:numel(f), f_stop, 'nearest' );
+f = f(f_idx_start:f_idx_stop);
+u = u(f_idx_start:f_idx_stop);
+i = i(f_idx_start:f_idx_stop);
+
+Z = abs(u./i);
+
+% analytic formular for characteristic impedance
+Z0 = sqrt(MUE0/EPS0) * log(coax_rad_ai/coax_rad_i) / (2*pi);
+upper_limit = Z0 * (1+upper_error);
+lower_limit = Z0 * (1-lower_error);
+
+if ENABLE_PLOTS
+ upper = upper_limit * ones(1,size(Z,2));
+ lower = lower_limit * ones(1,size(Z,2));
+ Z0_plot = Z0 * ones(1,size(Z,2));
+ figure
+ plot(f/1e9,[Z;upper;lower])
+ hold on
+ plot(f/1e9,Z0_plot,'m-.','LineWidth',2)
+ hold off
+ xlabel('Frequency (GHz)')
+ ylabel('Impedance (Ohm)')
+ legend( {'sim', 'upper limit', 'lower limit', 'theoretical'} );
+end
+
+pass = check_limits( Z, upper_limit, lower_limit );
+if pass
+ disp( 'combinedtests/Coax.m (characteristic impedance): pass' );
+else
+ disp( 'combinedtests/Coax.m (characteristic impedance): * FAILED *' );
+end
+
+
+
+
+if pass && CLEANUP
+ rmdir( [Sim_Path '/' Sim_CSX], 's' );
+end
+if ~pass && STOP_IF_FAILED
+ error 'test failed';
+end
diff --git a/TESTSUITE/combinedtests/README b/TESTSUITE/combinedtests/README
new file mode 100644
index 0000000..cc29c5b
--- /dev/null
+++ b/TESTSUITE/combinedtests/README
@@ -0,0 +1,3 @@
+#
+# These scripts test the full simulator (not single features)
+#
\ No newline at end of file
diff --git a/TESTSUITE/combinedtests/cavity.m b/TESTSUITE/combinedtests/cavity.m
new file mode 100644
index 0000000..39cbae9
--- /dev/null
+++ b/TESTSUITE/combinedtests/cavity.m
@@ -0,0 +1,212 @@
+function pass = cavity
+
+physical_constants;
+
+
+ENABLE_PLOTS = 1;
+CLEANUP = 0; % if enabled and result is PASS, remove simulation folder
+STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
+
+% LIMITS - inside
+lower_rel_limit = 1.3e-3; % -0.13%
+upper_rel_limit = 1.3e-3; % +0.13%
+lower_rel_limit_TM = 2.5e-3; % -0.25%
+upper_rel_limit_TM = 0; % +0%
+min_rel_amplitude = 0.6; % 60%
+min_rel_amplitude_TM = 0.27; % 27%
+
+% LIMITS - outside
+outer_rel_limit = 0.02;
+max_rel_amplitude = 0.17;
+
+
+% structure
+a = 5e-2;
+b = 2e-2;
+d = 6e-2;
+if ~((b f2, f2_idx = f2_idx - 1; end
+
+ if strcmp( type, 'inside' )
+ if max( val(f1_idx:f2_idx) ) < max1 * rel_amplitude
+ pass = false;
+ return
+ end
+ elseif strcmp( type, 'outside' )
+ if max( val(f1_idx:f2_idx) ) > max1 * rel_amplitude
+ pass = false;
+ return
+ end
+ else
+ error 'unsupported operation'
+ end
+end
diff --git a/TESTSUITE/helperscripts/check_limits.m b/TESTSUITE/helperscripts/check_limits.m
new file mode 100644
index 0000000..c1cba3d
--- /dev/null
+++ b/TESTSUITE/helperscripts/check_limits.m
@@ -0,0 +1,22 @@
+function pass = check_limits( Z, upper_limit, lower_limit )
+
+% make row vector
+if size(Z,1) ~= 1
+ Z = Z.';
+end
+
+if numel(upper_limit) == 1
+ upper_limit = upper_limit * ones(1,size(Z,2));
+end
+if numel(lower_limit) == 1
+ lower_limit = lower_limit * ones(1,size(Z,2));
+end
+
+
+pass = 1;
+if any( Z > upper_limit )
+ pass = 0;
+end
+if any( Z < lower_limit )
+ pass = 0;
+end
diff --git a/TESTSUITE/helperscripts/invoke_openEMS.m b/TESTSUITE/helperscripts/invoke_openEMS.m
new file mode 100644
index 0000000..e1a4c28
--- /dev/null
+++ b/TESTSUITE/helperscripts/invoke_openEMS.m
@@ -0,0 +1,15 @@
+function invoke_openEMS( opts )
+
+if nargin < 1
+ opts = '';
+end
+% openEMS_opts = [openEMS_opts ' --disable-dumps'];
+% openEMS_opts = [openEMS_opts ' --debug-material'];
+
+filename = mfilename('fullpath');
+dir = fileparts( filename );
+openEMS_Path = [dir '/../../'];
+
+command = [openEMS_Path 'openEMS.sh ' opts];
+disp(command);
+system(command);
diff --git a/TESTSUITE/helperscripts/physical_constants.m b/TESTSUITE/helperscripts/physical_constants.m
new file mode 100644
index 0000000..082efda
--- /dev/null
+++ b/TESTSUITE/helperscripts/physical_constants.m
@@ -0,0 +1,8 @@
+%
+% physical constants
+%
+
+% Bronstein 3rd ed., 1997, pp. 945-946
+c0 = 299792458; % m/s
+MUE0 = 4e-7*pi; % N/A^2
+EPS0 = 1/(MUE0*c0^2); % F/m
diff --git a/matlab/SetBoundaryCond.m b/matlab/SetBoundaryCond.m
index 6d1ba6d..ed5ba8c 100644
--- a/matlab/SetBoundaryCond.m
+++ b/matlab/SetBoundaryCond.m
@@ -1,4 +1,8 @@
function FDTD = SetBoundaryCond(FDTD,BC)
+% FDTD = SetBoundaryCond(FDTD,BC)
+%
+% BC = [xmin xmax ymin ymax zmin zmax];
+% ?min/?max: 0=PEC 1=PMC
FDTD.BoundaryCond.ATTRIBUTE.xmin=BC(1);
FDTD.BoundaryCond.ATTRIBUTE.xmax=BC(2);
diff --git a/openEMS.pro b/openEMS.pro
index 43c0d75..34c098e 100644
--- a/openEMS.pro
+++ b/openEMS.pro
@@ -1,7 +1,7 @@
# -------------------------------------------------
# Project created by QtCreator 2010-02-26T22:34:51
# -------------------------------------------------
-QT -= gui
+QT -= gui core
TARGET = openEMS
CONFIG += console
CONFIG -= app_bundle
@@ -14,11 +14,11 @@ LIBS += -L../CSXCAD \
-L../fparser \
-lfparser \
-L../tinyxml \
- -ltinyxml
-QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../CSXCAD\'
-QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../fparser\'
-QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../tinyxml\'
-
+ -ltinyxml \
+ -lboost_thread
+QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../CSXCAD\'
+QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../fparser\'
+QMAKE_LFLAGS += \'-Wl,-rpath,\$$ORIGIN/../tinyxml\'
SOURCES += main.cpp \
tools/ErrorMsg.cpp \
tools/AdrOp.cpp \
@@ -31,7 +31,8 @@ SOURCES += main.cpp \
FDTD/processfields_td.cpp \
FDTD/processcurrent.cpp \
examples/FDTD_examples.cpp \
- openems.cpp
+ openems.cpp \
+ FDTD/engine_multithread.cpp
HEADERS += tools/ErrorMsg.h \
tools/AdrOp.h \
tools/constants.h \
@@ -44,4 +45,8 @@ HEADERS += tools/ErrorMsg.h \
FDTD/processfields_td.h \
FDTD/processcurrent.h \
examples/FDTD_examples.h \
- openems.h
+ openems.h \
+ FDTD/engine_multithread.h
+
+QMAKE_CXXFLAGS_RELEASE = -O2 -g -march=native
+QMAKE_CXXFLAGS_DEBUG = -O0 -g -march=native
diff --git a/openems.cpp b/openems.cpp
index 1c311ab..2beb857 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -20,6 +20,7 @@
#include "tools/array_ops.h"
#include "FDTD/operator.h"
#include "FDTD/engine.h"
+#include "FDTD/engine_multithread.h"
#include "FDTD/processvoltage.h"
#include "FDTD/processcurrent.h"
#include "FDTD/processfields_td.h"
@@ -32,8 +33,8 @@
double CalcDiffTime(timeval t1, timeval t2)
{
- double s_diff = difftime(t1.tv_sec,t2.tv_sec);
- s_diff += ((double)t1.tv_usec-(double)t2.tv_usec)*1e-6;
+ double s_diff = t1.tv_sec - t2.tv_sec;
+ s_diff += (t1.tv_usec-t2.tv_usec)*1e-6;
return s_diff;
}
@@ -46,27 +47,22 @@ openEMS::openEMS()
DebugMat = false;
DebugOp = false;
endCrit = 1e-6;
+
+ m_engine = EngineType_Standard;
+ m_engine_numThreads = 0;
}
openEMS::~openEMS()
{
- delete FDTD_Eng;
- FDTD_Eng=NULL;
- delete PA;
- PA=NULL;
- delete FDTD_Op;
- FDTD_Op=NULL;
+ Reset();
}
void openEMS::Reset()
{
- delete FDTD_Op;
- FDTD_Op=NULL;
- delete FDTD_Eng;
- FDTD_Eng=NULL;
if (PA) PA->DeleteAll();
- delete PA;
- PA=NULL;
+ delete PA; PA=0;
+ delete FDTD_Eng; FDTD_Eng=0;
+ delete FDTD_Op; FDTD_Op=0;
}
//! \brief processes a command line argument
@@ -95,6 +91,18 @@ bool openEMS::parseCommandLineArgument( const char *argv )
DebugOperator();
return true;
}
+ else if (strcmp(argv,"--engine=multithreaded")==0)
+ {
+ cout << "openEMS - enabled multithreading" << endl;
+ m_engine = EngineType_Multithreaded;
+ return true;
+ }
+ else if (strncmp(argv,"--numThreads=",13)==0)
+ {
+ m_engine_numThreads = atoi(argv+13);
+ cout << "openEMS - fixed number of threads: " << m_engine_numThreads << endl;
+ return true;
+ }
return false;
}
@@ -158,6 +166,11 @@ int openEMS::SetupFDTD(const char* file)
Excite->QueryDoubleAttribute("f0",&f0);
fc = 0;
}
+ else if (Excit_Type==2)
+ {
+ Excite->QueryDoubleAttribute("f0",&f0);
+ fc = 0;
+ }
TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond");
if (BC==NULL)
@@ -210,6 +223,24 @@ int openEMS::SetupFDTD(const char* file)
exit(2);
}
}
+ else if (Excit_Type==2)
+ {
+ Nyquist = FDTD_Op->CalcDiracPulsExcitation();
+ if (!Nyquist)
+ {
+ cerr << "openEMS: excitation setup failed!!" << endl;
+ exit(2);
+ }
+ }
+ else if (Excit_Type==3)
+ {
+ Nyquist = FDTD_Op->CalcStepExcitation();
+ if (!Nyquist)
+ {
+ cerr << "openEMS: excitation setup failed!!" << endl;
+ exit(2);
+ }
+ }
else
{
cerr << "openEMS: Excitation type is unknown" << endl;
@@ -234,7 +265,15 @@ int openEMS::SetupFDTD(const char* file)
cout << "Creation time for operator: " << difftime(OpDoneTime,startTime) << " s" << endl;
//create FDTD engine
- FDTD_Eng = new Engine(FDTD_Op);
+ switch (m_engine) {
+ case EngineType_Multithreaded:
+ FDTD_Eng = Engine_Multithread::createEngine(FDTD_Op,m_engine_numThreads);
+ break;
+ default:
+ FDTD_Eng = Engine::createEngine(FDTD_Op);
+ break;
+ }
+
time_t currTime = time(NULL);
@@ -317,17 +356,20 @@ void openEMS::RunFDTD()
{
cout << "Running FDTD engine... this may take a while... grab a cup of coffee?!?" << endl;
- timeval currTime;
- gettimeofday(&currTime,NULL);
- timeval startTime = currTime;
- timeval prevTime= currTime;
ProcessFields ProcField(FDTD_Op,FDTD_Eng);
double maxE=0,currE=0;
double change=1;
int prevTS=0,currTS=0;
- double speed = (double)FDTD_Op->GetNumberCells()/1e6;
+ double speed = FDTD_Op->GetNumberCells()/1e6;
double t_diff;
+
+ timeval currTime;
+ gettimeofday(&currTime,NULL);
+ timeval startTime = currTime;
+ timeval prevTime= currTime;
+
//*************** simulate ************//
+
int step=PA->Process();
if ((step<0) || (step>(int)NrTS)) step=NrTS;
while ((FDTD_Eng->GetNumberOfTimesteps()endCrit))
@@ -347,7 +389,7 @@ void openEMS::RunFDTD()
if (currE>maxE)
maxE=currE;
cout << "[@" << setw(8) << (int)CalcDiffTime(currTime,startTime) << "s] Timestep: " << setw(12) << currTS << " (" << setw(6) << setprecision(2) << std::fixed << (double)currTS/(double)NrTS*100.0 << "%)" ;
- cout << " with currently " << setw(6) << setprecision(1) << std::fixed << speed*(double)(currTS-prevTS)/t_diff << " MCells/s" ;
+ cout << " with currently " << setw(6) << setprecision(1) << std::fixed << speed*(currTS-prevTS)/t_diff << " MCells/s" ;
if (maxE)
change = currE/maxE;
cout << " --- Energy: ~" << setw(6) << setprecision(2) << std::scientific << currE << " (decrement: " << setw(6) << setprecision(2) << std::fixed << fabs(10.0*log10(change)) << "dB)" << endl;
diff --git a/openems.h b/openems.h
index ca3378b..981ac30 100644
--- a/openems.h
+++ b/openems.h
@@ -52,6 +52,10 @@ protected:
Operator* FDTD_Op;
Engine* FDTD_Eng;
ProcessingArray* PA;
+
+ enum EngineType {EngineType_Standard,EngineType_Multithreaded};
+ EngineType m_engine;
+ unsigned int m_engine_numThreads;
};
#endif // OPENEMS_H
diff --git a/tools/array_ops.cpp b/tools/array_ops.cpp
index 42c010d..d80a231 100644
--- a/tools/array_ops.cpp
+++ b/tools/array_ops.cpp
@@ -18,7 +18,7 @@
#include "array_ops.h"
#include
-FDTD_FLOAT*** Create3DArray(unsigned int* numLines)
+FDTD_FLOAT*** Create3DArray(const unsigned int* numLines)
{
FDTD_FLOAT*** array=NULL;
unsigned int pos[3];
@@ -38,7 +38,7 @@ FDTD_FLOAT*** Create3DArray(unsigned int* numLines)
return array;
}
-void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines)
+void Delete3DArray(FDTD_FLOAT*** array, const unsigned int* numLines)
{
if (array==NULL) return;
unsigned int pos[3];
@@ -53,7 +53,7 @@ void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines)
delete[] array;
}
-FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines)
+FDTD_FLOAT**** Create_N_3DArray(const unsigned int* numLines)
{
FDTD_FLOAT**** array=NULL;
array = new FDTD_FLOAT***[3];
@@ -64,7 +64,7 @@ FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines)
return array;
}
-void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines)
+void Delete_N_3DArray(FDTD_FLOAT**** array, const unsigned int* numLines)
{
if (array==NULL) return;
unsigned int pos[3];
diff --git a/tools/array_ops.h b/tools/array_ops.h
index b04c821..5424e05 100644
--- a/tools/array_ops.h
+++ b/tools/array_ops.h
@@ -20,12 +20,12 @@
#include "../FDTD/operator.h"
-FDTD_FLOAT*** Create3DArray(unsigned int* numLines);
-void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines);
+FDTD_FLOAT*** Create3DArray(const unsigned int* numLines);
+void Delete3DArray(FDTD_FLOAT*** array, const unsigned int* numLines);
-FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines);
-void Delete_N_3DArray(FDTD_FLOAT**** array, unsigned int* numLines);
+FDTD_FLOAT**** Create_N_3DArray(const unsigned int* numLines);
+void Delete_N_3DArray(FDTD_FLOAT**** array, const unsigned int* numLines);
-void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, unsigned int* numLines);
+void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, const unsigned int* numLines);
#endif // ARRAY_OPS_H