first working MPI implementation

todo:
 - only/mainly rank 0 takes care of output
 - rank 0 manages end-criteria, energy calculation
 - field dumps split between ranks --> prevent file corruption
 - a lot more...

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

FDTD/engine_mpi.cpp

@@ -39,9 +39,165 @@ void Engine_MPI::Init()
 {
 	Engine_SSE_Compressed::Init();
 
+	for (int i=0;i<3;++i)
+	{
+		m_BufferUp[i]=NULL;
+		m_BufferDown[i]=NULL;
+		m_BufferSize[i]=0;
+	}
+
+	if (m_Op_MPI->m_NumProc>1)
+	{
+		// init buffers, nx*ny*2 for the tangential electric or magnetic fields at the interface
+		m_BufferSize[2] = m_Op_MPI->numLines[0]*m_Op_MPI->numLines[1]*2;
+		if (m_Op_MPI->m_NeighborDown[2]>=0)
+		{
+			m_BufferDown[2] = new float[m_BufferSize[2]];
+		}
+		if (m_Op_MPI->m_NeighborUp[2]>=0)
+		{
+			m_BufferUp[2] = new float[m_BufferSize[2]];
+		}
+	}
+
 }
 
 void Engine_MPI::Reset()
 {
+	for (int i=0;i<3;++i)
+	{
+		delete[] m_BufferUp[i];
+		delete[] m_BufferDown[i];
+		m_BufferUp[i]=NULL;
+		m_BufferDown[i]=NULL;
+		m_BufferSize[i]=0;
+	}
+
 	Engine_SSE_Compressed::Reset();
 }
+
+void Engine_MPI::SendReceiveVoltages()
+{
+	unsigned int pos[3];
+
+	//non-blocking prepare for receive...
+	for (int n=0;n<3;++n)
+		if (m_Op_MPI->m_NeighborDown[n]>=0)
+			MPI_Irecv( m_BufferDown[n] , m_BufferSize[n], MPI_FLOAT, m_Op_MPI->m_NeighborDown[n], m_Op_MPI->m_MyTag, MPI_COMM_WORLD, &Recv_Request[n]);
+
+	for (int n=0;n<3;++n)
+	{
+		int nP  = (n+1)%3;
+		int nPP = (n+2)%3;
+
+		//send voltages
+		unsigned int iPos=0;
+		pos[n]=numLines[n]-1;
+		if (m_Op_MPI->m_NeighborUp[n]>=0)
+		{
+			for (pos[nP]=0; pos[nP]<numLines[nP]-1; ++pos[nP])
+			{
+				for (pos[nPP]=0; pos[nPP]<numLines[nPP]-1; ++pos[nPP])
+				{
+					m_BufferUp[n][iPos++] = Engine_SSE_Compressed::GetVolt(nP ,pos);
+					m_BufferUp[n][iPos++] = Engine_SSE_Compressed::GetVolt(nPP,pos);
+				}
+			}
+			MPI_Isend( m_BufferUp[n] , m_BufferSize[n], MPI_FLOAT, m_Op_MPI->m_NeighborUp[n], m_Op_MPI->m_MyTag, MPI_COMM_WORLD, &Send_Request[n]);
+		}
+
+		//receive voltages
+		pos[n]=0;
+		iPos=0;
+		if (m_Op_MPI->m_NeighborDown[n]>=0)
+		{
+			//wait for receive to finish...
+			MPI_Wait(&Recv_Request[n],&stat);
+			for (pos[nP]=0; pos[nP]<numLines[nP]-1; ++pos[nP])
+			{
+				for (pos[nPP]=0; pos[nPP]<numLines[nPP]-1; ++pos[nPP])
+				{
+					Engine_SSE_Compressed::SetVolt(nP ,pos,m_BufferDown[n][iPos++]);
+					Engine_SSE_Compressed::SetVolt(nPP,pos,m_BufferDown[n][iPos++]);
+				}
+			}
+		}
+
+	}
+}
+
+void Engine_MPI::SendReceiveCurrents()
+{
+	unsigned int pos[3];
+
+	//non-blocking prepare for receive...
+	for (int n=0;n<3;++n)
+		if (m_Op_MPI->m_NeighborUp[n]>=0)
+			MPI_Irecv( m_BufferUp[n] , m_BufferSize[n], MPI_FLOAT, m_Op_MPI->m_NeighborUp[n], m_Op_MPI->m_MyTag, MPI_COMM_WORLD, &Recv_Request[n]);
+
+	for (int n=0;n<3;++n)
+	{
+		int nP  = (n+1)%3;
+		int nPP = (n+2)%3;
+
+		//send currents
+		unsigned int iPos=0;
+		pos[n]=0;
+		if (m_Op_MPI->m_NeighborDown[n]>=0)
+		{
+			for (pos[nP]=0; pos[nP]<numLines[nP]-1; ++pos[nP])
+			{
+				for (pos[nPP]=0; pos[nPP]<numLines[nPP]-1; ++pos[nPP])
+				{
+					m_BufferDown[n][iPos++] = Engine_SSE_Compressed::GetCurr(nP ,pos);
+					m_BufferDown[n][iPos++] = Engine_SSE_Compressed::GetCurr(nPP,pos);
+				}
+			}
+			MPI_Isend( m_BufferDown[n] , m_BufferSize[n], MPI_FLOAT, m_Op_MPI->m_NeighborDown[n], m_Op_MPI->m_MyTag, MPI_COMM_WORLD, &Send_Request[n]);
+		}
+
+		//receive currents
+		pos[n]=numLines[n]-1;
+		iPos=0;
+		if (m_Op_MPI->m_NeighborUp[n]>=0)
+		{
+			//wait for receive to finish...
+			MPI_Wait(&Recv_Request[n],&stat);
+			for (pos[nP]=0; pos[nP]<numLines[nP]-1; ++pos[nP])
+			{
+				for (pos[nPP]=0; pos[nPP]<numLines[nPP]-1; ++pos[nPP])
+				{
+					Engine_SSE_Compressed::SetCurr(nP ,pos,m_BufferUp[n][iPos++]);
+					Engine_SSE_Compressed::SetCurr(nPP,pos,m_BufferUp[n][iPos++]);
+				}
+			}
+		}
+
+	}
+}
+
+bool Engine_MPI::IterateTS(unsigned int iterTS)
+{
+	for (unsigned int iter=0; iter<iterTS; ++iter)
+	{
+		//voltage updates with extensions
+		DoPreVoltageUpdates();
+		UpdateVoltages(0,numLines[0]);
+		DoPostVoltageUpdates();
+		Apply2Voltages();
+		ApplyVoltageExcite();
+		SendReceiveVoltages();
+
+		//current updates with extensions
+		DoPreCurrentUpdates();
+		UpdateCurrents(0,numLines[0]-1);
+		DoPostCurrentUpdates();
+		Apply2Current();
+		ApplyCurrentExcite();
+		SendReceiveCurrents();
+
+		++numTS;
+	}
+	return true;
+}
+

FDTD/engine_mpi.h

@@ -20,6 +20,7 @@
 
 #include "operator_mpi.h"
 #include "engine_sse_compressed.h"
+#include "mpi.h"
 
 class Engine_MPI : public Engine_SSE_Compressed
 {
@@ -32,10 +33,23 @@ public:
 	virtual void Init();
 	virtual void Reset();
 
+	virtual bool IterateTS(unsigned int iterTS);
+
 protected:
 	Engine_MPI(const Operator_MPI* op);
 	const Operator_MPI* m_Op_MPI;
 
+	MPI_Status stat;
+	MPI_Request Send_Request[3];
+	MPI_Request Recv_Request[3];
+
+	//field buffer for MPI transfer...
+	unsigned int m_BufferSize[3];
+	float* m_BufferUp[3];
+	float* m_BufferDown[3];
+
+	void SendReceiveVoltages();
+	void SendReceiveCurrents();
 };
 
 #endif // ENGINE_MPI_H

FDTD/operator_mpi.cpp

@@ -16,11 +16,12 @@
 */
 
 #include "operator_mpi.h"
-
 #include "operator_sse_compressed.h"
 #include "engine_sse_compressed.h"
 #include "engine_mpi.h"
 #include "tools/array_ops.h"
+#include "tools/useful.h"
+#include "mpi.h"
 
 Operator_MPI* Operator_MPI::New()
 {
@@ -39,23 +40,98 @@ Operator_MPI::~Operator_MPI()
 	Delete();
 }
 
+bool Operator_MPI::SetGeometryCSX(ContinuousStructure* geo)
+{
+	//manipulate geometry for this part...
+
+	if (m_NumProc>1)
+	{
+		CSRectGrid* grid = geo->GetGrid();
+		int nz = grid->GetQtyLines(2);
+		std::vector<unsigned int> jobs = AssignJobs2Threads(nz, m_NumProc);
+		double z_lines[jobs.at(m_MyID)];
+		unsigned int z_start=0;
+		for (int n=0;n<m_MyID;++n)
+			z_start+=jobs.at(n);
+		for (unsigned int n=0;n<jobs.at(m_MyID);++n)
+			z_lines[n] = grid->GetLine(2,z_start+n);
+		grid->ClearLines(2);
+		grid->AddDiscLines(2,jobs.at(m_MyID),z_lines);
+
+		//lower neighbor is ID-1
+		if (m_MyID>0)
+			m_NeighborDown[2]=m_MyID-1;
+		//upper neighbor is ID+1
+		if (m_MyID<m_NumProc-1)
+			m_NeighborUp[2]=m_MyID+1;
+	}
+	else
+		cerr << "Operator_MPI::SetGeometryCSX: Warning: Number of MPI processes is 1, skipping MPI engine... " << endl;
+
+	return Operator_SSE_Compressed::SetGeometryCSX(geo);
+}
+
+void Operator_MPI::SetBoundaryCondition(int* BCs)
+{
+	//set boundary conditions on MPI interfaces to PEC, ApplyElectricBC will handle proper interface handling...
+	for (int n=0;n<3;++n)
+	{
+		if (m_NeighborUp[n]>=0)
+			BCs[2*n+1] = 0;
+		if (m_NeighborDown[n]>=0)
+			BCs[2*n] = 0;
+	}
+	Operator_SSE_Compressed::SetBoundaryCondition(BCs);
+}
+
+void Operator_MPI::ApplyElectricBC(bool* dirs)
+{
+	for (int n=0;n<3;++n)
+	{
+		//do not delete operator at upper inteface
+		if (m_NeighborUp[n]>=0)
+			dirs[2*n+1] = false;
+	}
+	Operator_SSE_Compressed::ApplyElectricBC(dirs);
+}
+
 Engine* Operator_MPI::CreateEngine() const
 {
-	Engine_MPI* eng = Engine_MPI::New(this);
+	if (m_NumProc>1)
-	return eng;
+		return Engine_MPI::New(this);
+	else
+		return Engine_SSE_Compressed::New(this);
 }
 
 void Operator_MPI::Init()
 {
-	Operator_sse::Init();
+	Operator_SSE_Compressed::Init();
+
+	m_MyTag = 0;
+
+	for (int i=0;i<3;++i)
+	{
+		m_NeighborUp[i]=-1;
+		m_NeighborDown[i]=-1;
+	}
+
+	int  namelen;
+	m_NumProc = MPI::COMM_WORLD.Get_size();
+	m_MyID = MPI::COMM_WORLD.Get_rank();
+
+	m_Processor_Name = new char[MPI_MAX_PROCESSOR_NAME];
+	MPI::Get_processor_name(m_Processor_Name,namelen);
+
+	cerr << "Operator_MPI::Init(): Running on " << m_Processor_Name << endl;
 }
 
 void Operator_MPI::Delete()
 {
+	delete[] m_Processor_Name;
 }
 
 void Operator_MPI::Reset()
 {
 	Delete();
-	Operator_sse::Reset();
+	Operator_SSE_Compressed::Reset();
 }

FDTD/operator_mpi.h

@@ -22,11 +22,17 @@
 
 class Operator_MPI : public Operator_SSE_Compressed
 {
+	friend class Engine_MPI;
 public:
 	//! Create a new operator
 	static Operator_MPI* New();
 	virtual ~Operator_MPI();
 
+	virtual bool SetGeometryCSX(ContinuousStructure* geo);
+
+	virtual void SetBoundaryCondition(int* BCs);
+	virtual void ApplyElectricBC(bool* dirs);
+
 	virtual Engine* CreateEngine() const;
 
 protected:
@@ -34,6 +40,15 @@ protected:
 	virtual void Init();
 	void Delete();
 	virtual void Reset();
+
+	int m_MyID;
+	int m_NumProc;
+	int m_MyTag;
+	char* m_Processor_Name;
+
+	//the up and down neighbors, -1 if non for the given direction
+	int m_NeighborUp[3];
+	int m_NeighborDown[3];
 };
 
 #endif // OPERATOR_MPI_H