Merge branch 'multithreading'

Adding multithreading support to openEMS

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

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;

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 <http://www.gnu.org/licenses/>.
+*/
+
+//#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 <iomanip>
+
+//! \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<boost::thread::id, std::vector<double> >::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<double>::iterator it2;
+		for (it2=it->second.begin(); it2<it->second.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; n<m_numThreads; n++) {
+		unsigned int start = n * linesPerThread;
+		unsigned int stop = (n+1) * linesPerThread - 1;
+		unsigned int stop_h = stop;
+		if (n == m_numThreads-1) {
+			// last thread
+			stop = Op->numLines[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 <iterTS> 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;iter<m_enginePtr->m_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]<m_enginePtr->Op->numLines[1];++pos[1])
+				{
+					shift[1]=pos[1];
+					for (pos[2]=0;pos[2]<m_enginePtr->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]<m_enginePtr->Op->numLines[1]-1;++pos[1])
+				{
+					for (pos[2]=0;pos[2]<m_enginePtr->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;n<E_Exc_Count;++n)
+		{
+			exc_pos = (int)m_enginePtr->numTS - (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

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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef ENGINE_MULTITHREAD_H
+#define ENGINE_MULTITHREAD_H
+
+#include "operator.h"
+#include "engine.h"
+
+#include <boost/thread.hpp>
+#include <boost/fusion/include/list.hpp>
+#include <boost/fusion/container/list/list_fwd.hpp>
+#include <boost/fusion/include/list_fwd.hpp>
+
+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<boost::thread::id, std::vector<double> > m_timer_list;
+#endif
+};
+
+#endif // ENGINE_MULTITHREAD_H

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;

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

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

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

TESTSUITE/combinedtests/README

@@ -0,0 +1,3 @@
+#
+# These scripts test the full simulator (not single features)
+#

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<a) && (a<d))
+    error 'correct the dimensions of the cavity'
+end
+
+f_start = 1e9;
+f_stop = 10e9;
+
+Sim_Path = 'tmp';
+Sim_CSX = 'cavity.xml';
+
+[status,message,messageid]=mkdir(Sim_Path);
+
+%setup FDTD parameter
+FDTD = InitFDTD( 10000,1e-6 );
+FDTD = SetGaussExcite(FDTD,(f_stop-f_start)/2,(f_stop-f_start)/2);
+BC = [0 0 0 0 0 0]; % PEC boundaries
+FDTD = SetBoundaryCond(FDTD,BC);
+
+%setup CSXCAD geometry
+CSX = InitCSX();
+grid_res = 2e-3;
+mesh.x = 0:grid_res:a; %linspace(0,a,25);
+mesh.y = 0:grid_res:b; %linspace(0,b,25);
+mesh.z = 0:grid_res:d; %linspace(0,d,25);
+CSX = DefineRectGrid(CSX, 1,mesh);
+
+% excitation
+CSX = AddExcitation(CSX,'excite1',0,[1 1 1]);
+p(1,1) = mesh.x(floor(end*2/3));
+p(2,1) = mesh.y(floor(end*2/3));
+p(3,1) = mesh.z(floor(end*2/3));
+p(1,2) = mesh.x(floor(end*2/3)+1);
+p(2,2) = mesh.y(floor(end*2/3)+1);
+p(3,2) = mesh.z(floor(end*2/3)+1);
+CSX = AddCurve( CSX, 'excite1', 0, p );
+ 
+%dump
+% CSX = AddDump(CSX,'Et_',0,2);
+% pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(end) mesh.z(end)];
+% CSX = AddBox(CSX,'Et_',0 , pos1,pos2);
+
+% %dump
+% CSX = AddDump(CSX,'Et2_',0,2);
+% pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(1) mesh.z(end)];
+% CSX = AddBox(CSX,'Et2_',0 , pos1,pos2);
+% 
+% %dump
+% CSX = AddDump(CSX,'Et3_',0,2);
+% pos1 = [mesh.x(1) mesh.y(end-1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(end-1) mesh.z(end)];
+% CSX = AddBox(CSX,'Et3_',0 , pos1,pos2);
+
+%voltage calc
+CSX = AddProbe(CSX,'ut1x',0);
+pos1 = [mesh.x(floor(end/4)) mesh.y(floor(end/2))   mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/4)+1) mesh.y(floor(end/2)) mesh.z(floor(end/5))];
+CSX = AddBox(CSX,'ut1x', 0 ,pos1,pos2);
+
+CSX = AddProbe(CSX,'ut1y',0);
+pos1 = [mesh.x(floor(end/4)) mesh.y(floor(end/2))   mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/4)) mesh.y(floor(end/2)+1) mesh.z(floor(end/5))];
+CSX = AddBox(CSX,'ut1y', 0 ,pos1,pos2);
+
+CSX = AddProbe(CSX,'ut1z',0);
+pos1 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(floor(end/5)+1)];
+CSX = AddBox(CSX,'ut1z', 0 ,pos1,pos2);
+
+%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( {'ut1x', 'ut1y', 'ut1z'} );
+cd(savePath);
+
+
+
+%
+% analysis
+%
+
+% remove excitation from time series
+t_start = 7e-10; % FIXME to be calculated
+t_idx_start = interp1( UI.TD{1}.t, 1:numel(UI.TD{1}.t), t_start, 'nearest' );
+for n=1:numel(UI.TD)
+    UI.TD{n}.t = UI.TD{n}.t(t_idx_start:end);
+    UI.TD{n}.val = UI.TD{n}.val(t_idx_start:end);
+    [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( UI.TD{n}.t, UI.TD{n}.val );
+end
+
+
+f = UI.FD{1}.f;
+ux = UI.FD{1}.val;
+uy = UI.FD{2}.val;
+uz = UI.FD{3}.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);
+ux = ux(f_idx_start:f_idx_stop);
+uy = uy(f_idx_start:f_idx_stop);
+uz = uz(f_idx_start:f_idx_stop);
+
+% analytic formula for resonant wavenumber
+k = @(m,n,l) sqrt( (m*pi/a)^2 + (n*pi/b)^2 + (l*pi/d)^2 );
+f_TE101 = c0/(2*pi) * k(1,0,1);
+f_TE102 = c0/(2*pi) * k(1,0,2);
+f_TE201 = c0/(2*pi) * k(2,0,1);
+f_TE202 = c0/(2*pi) * k(2,0,2);
+f_TM110 = c0/(2*pi) * k(1,1,0);
+f_TM111 = c0/(2*pi) * k(1,1,1);
+
+f_TE = [f_TE101 f_TE102 f_TE201 f_TE202];
+f_TM = [f_TM110 f_TM111];
+
+% calculate frequency limits
+temp = [f_start f_TE f_stop];
+f_outer1 = [];
+f_outer2 = [];
+for n=1:numel(temp)-1
+    f_outer1 = [f_outer1 temp(n) .* (1+outer_rel_limit)];
+    f_outer2 = [f_outer2 temp(n+1) .* (1-outer_rel_limit)];
+end
+
+temp = [f_start f_TM f_stop];
+f_outer1_TM = [];
+f_outer2_TM = [];
+for n=1:numel(temp)-1
+    f_outer1_TM = [f_outer1_TM temp(n) .* (1+outer_rel_limit)];
+    f_outer2_TM = [f_outer2_TM temp(n+1) .* (1-outer_rel_limit)];
+end
+
+
+if ENABLE_PLOTS
+    figure
+    plot(f/1e9,abs(uy))
+    max1 = max(abs(uy));
+    hold on
+    plot( repmat(f_TE,2,1)/1e9, repmat([0; max1],1,numel(f_TE)), 'm-.', 'LineWidth', 2 )
+    plot( (repmat(f_TE,2,1) .* repmat(1-lower_rel_limit,2,numel(f_TE)))/1e9, repmat([0; max1],1,numel(f_TE)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TE,2,1) .* repmat(1+upper_rel_limit,2,numel(f_TE)))/1e9, repmat([0; max1],1,numel(f_TE)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TE,2,1) .* repmat([1-outer_rel_limit;1+outer_rel_limit],1,numel(f_TE)))/1e9, repmat(max1*min_rel_amplitude,2,numel(f_TE)), 'r-', 'LineWidth', 1 ) % freq limits
+    plot( [f_outer1;f_outer2]/1e9, repmat(max1*max_rel_amplitude,2,numel(f_outer1)), 'g-', 'LineWidth', 1 ) % amplitude limits
+    xlabel('Frequency (GHz)')
+    legend( {'u_y','theoretical'} )
+    title( 'TE-modes' )
+
+    figure
+    plot(f/1e9,abs(uz))
+    max1 = max(abs(uz));
+    hold on
+    plot( repmat(f_TM,2,1)/1e9, repmat([0; max1],1,numel(f_TM)), 'm-.', 'LineWidth', 2 )
+    plot( (repmat(f_TM,2,1) .* repmat(1-lower_rel_limit_TM,2,numel(f_TM)))/1e9, repmat([0; max1],1,numel(f_TM)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TM,2,1) .* repmat(1+upper_rel_limit_TM,2,numel(f_TM)))/1e9, repmat([0; max1],1,numel(f_TM)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TM,2,1) .* repmat([1-lower_rel_limit_TM;1+upper_rel_limit_TM],1,numel(f_TM)))/1e9, repmat(max1*min_rel_amplitude_TM,2,numel(f_TM)), 'r-', 'LineWidth', 1 ) % freq limits
+    plot( [f_outer1_TM;f_outer2_TM]/1e9, repmat(max1*max_rel_amplitude,2,numel(f_outer1_TM)), 'g-', 'LineWidth', 1 ) % amplitude limits
+    xlabel('Frequency (GHz)')
+    legend( {'u_z','theoretical'} )
+    title( 'TM-modes' )
+end
+
+pass1 = check_frequency( f, abs(uy), f_TE*(1+upper_rel_limit),    f_TE*(1-lower_rel_limit),    min_rel_amplitude, 'inside' );
+pass2 = check_frequency( f, abs(uz), f_TM*(1+upper_rel_limit_TM), f_TM*(1-lower_rel_limit_TM), min_rel_amplitude_TM, 'inside' );
+pass3 = check_frequency( f, abs(uy), f_outer2, f_outer1, max_rel_amplitude, 'outside' );
+pass4 = check_frequency( f, abs(uz), f_outer2_TM, f_outer1_TM, max_rel_amplitude, 'outside' );
+pass = pass1 && pass2 && pass3 && pass4;
+if pass
+    disp( 'combinedtests/cavity.m (resonance frequency):  pass' );
+else
+    disp( 'combinedtests/cavity.m (resonance frequency):  * FAILED *' );
+end
+
+
+
+
+if pass && CLEANUP
+    rmdir( [Sim_Path '/' Sim_CSX], 's' );
+end
+if ~pass && STOP_IF_FAILED
+    error 'test failed';
+end

TESTSUITE/helperscripts/check_frequency.m

@@ -0,0 +1,31 @@
+function pass = check_frequency( f, val, f_upper, f_lower, rel_amplitude, type )
+
+pass = true;
+max1 = max(val);
+
+if numel(f_upper) ~= numel(f_lower)
+    error 'inconsistant vectors'
+end
+
+for n=1:numel(f_upper)
+    f1 = f_lower(n);
+    f2 = f_upper(n);
+    f1_idx = interp1( f, 1:numel(f), f1, 'nearest' );
+%    if f(f1_idx) < f1, f1_idx = f1_idx + 1; end
+    f2_idx = interp1( f, 1:numel(f), f2, 'nearest' );
+%    if f(f2_idx) > 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

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

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

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

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

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

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);
+	double s_diff = t1.tv_sec - t2.tv_sec;
-	s_diff += ((double)t1.tv_usec-(double)t2.tv_usec)*1e-6;
+	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;
+	Reset();
-	FDTD_Eng=NULL;
-	delete PA;
-	PA=NULL;
-	delete FDTD_Op;
-	FDTD_Op=NULL;
 }
 
 void openEMS::Reset()
 {
-	delete FDTD_Op;
-	FDTD_Op=NULL;
-	delete FDTD_Eng;
-	FDTD_Eng=NULL;
 	if (PA) PA->DeleteAll();
-	delete PA;
+	delete PA; PA=0;
-	PA=NULL;
+	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()<NrTS) && (change>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;

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

tools/array_ops.cpp

@@ -18,7 +18,7 @@
 #include "array_ops.h"
 #include <ostream>
 
-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];

tools/array_ops.h

@@ -20,12 +20,12 @@
 
 #include "../FDTD/operator.h"
 
-FDTD_FLOAT*** Create3DArray(unsigned int* numLines);
+FDTD_FLOAT*** Create3DArray(const unsigned int* numLines);
-void Delete3DArray(FDTD_FLOAT*** array, unsigned int* numLines);
+void Delete3DArray(FDTD_FLOAT*** array, const unsigned int* numLines);
 
-FDTD_FLOAT**** Create_N_3DArray(unsigned int* numLines);
+FDTD_FLOAT**** Create_N_3DArray(const unsigned int* numLines);
-void Delete_N_3DArray(FDTD_FLOAT**** array, 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