diff --git a/bin/beyn.py b/bin/beyn.py
index 9f9c2a573f0a7ddd29dda5058c24d74a495cfdf2..c0628d62c173f970c45878817d68e23fccde5f05 100644
--- a/bin/beyn.py
+++ b/bin/beyn.py
@@ -7,6 +7,7 @@ Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
+from mpi4py import MPI
import solver as Solver
import sys as sys
import numpy as np
@@ -42,7 +43,9 @@ class Beyn:
self.__lStart = lStart
self.__lStep = lStep
self.__rankTol = rankTol
- self.__verbose = verbose
+ self.__myProc = MPI.COMM_WORLD.Get_rank()
+ self.__nProc = MPI.COMM_WORLD.Get_size()
+ self.__verbose = verbose and (self.__myProc == 0)
def simple(self):
"""Solves the eigenvalue problem using Beyn's algorithm (simple version)
@@ -137,10 +140,14 @@ class Beyn:
If beQuite is True, this solver becomes quite; it is verbose otherwise
"""
- self.__verbose = not beQuite
+ if(self.__myProc == 0):
+ self.__verbose = not beQuite
def display(self):
"""Prints this Beyn solver parameters"""
+ if(self.__myProc != 0):
+ return
+
print "Beyn's contour integral method (simple)"
print "---------------------------------------"
print " # Nodes used for the trapezoidal rule:"+ " " +str(self.__nNode)
@@ -177,18 +184,38 @@ class Beyn:
vHat -- the RHS matrix defing Beyn's integrals
"""
- # Initialise integrals
- A0 = np.matlib.zeros(vHat.shape, dtype=complex)
- A1 = np.matlib.zeros(vHat.shape, dtype=complex)
+ # Share integration points among processes
+ baseShare = self.__nNode // self.__nProc #Base share for proc
+ residualShare = self.__nNode % self.__nProc #Residual share
+ allStart = [0] #Inital start
+
+ for i in range(residualShare): #Proc with extra share
+ allStart.append(allStart[-1] + baseShare + 1) # -> previous+base+1
+ for i in range(residualShare, self.__nProc): #Proc with base share
+ allStart.append(allStart[-1] + baseShare) # -> previous+base+0
+
+ myStart = allStart[self.__myProc] #This proc start
+ myStop = allStart[self.__myProc+1] #This proc end
+
+ # Initialise integrals (partial results)
+ A0_partial = np.matlib.zeros(vHat.shape, dtype=complex)
+ A1_partial = np.matlib.zeros(vHat.shape, dtype=complex)
# Loop over integation points and integrate
- for i in range(self.__nNode):
+ for i in range(myStart, myStop):
t = 2 * np.pi * i / self.__nNode;
phi = self.__origin + self.__radius * np.exp(1j * t)
tmp = self.__multiSolve(vHat, phi)
- A0 += tmp * np.power(phi, 0) * np.exp(1j * t)
- A1 += tmp * np.power(phi, 1) * np.exp(1j * t)
+ A0_partial += tmp * np.power(phi, 0) * np.exp(1j * t)
+ A1_partial += tmp * np.power(phi, 1) * np.exp(1j * t)
+ MPI.COMM_WORLD.Barrier()
+
+ # Synchronize parital results
+ A0 = np.matlib.zeros(vHat.shape, dtype=complex)
+ A1 = np.matlib.zeros(vHat.shape, dtype=complex)
+ MPI.COMM_WORLD.Allreduce(A0_partial, A0, op=MPI.SUM)
+ MPI.COMM_WORLD.Allreduce(A1_partial, A1, op=MPI.SUM)
# Final scale
A0 *= self.__radius / self.__nNode
@@ -204,7 +231,6 @@ class Beyn:
B -- the matrix of RHS
w -- the complex frequency to use
"""
-
# Solve
self.__operator.solve(B, w)
@@ -218,7 +244,13 @@ class Beyn:
# Done
return X
- @staticmethod
- def __randomMatrix(n, m):
+ def __randomMatrix(self, n, m):
"""Returns a random complex matrix of the given size (n x m)"""
- return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
+ if(self.__myProc == 0):
+ A = np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
+ else:
+ A = None
+
+ A = MPI.COMM_WORLD.bcast(A, root=0)
+ MPI.COMM_WORLD.Barrier()
+ return A
diff --git a/bin/cim.py b/bin/cim.py
index 2ab24cf3a21c0f7e5e5d8bbdf0134604680c1030..142a838829b59a7fd6e9e989255ddc56082d8781 100755
--- a/bin/cim.py
+++ b/bin/cim.py
@@ -11,6 +11,7 @@ See the LICENSE.txt and README.md for more license and copyright information.
from solver import GetDPWave
from beyn import Beyn
from fixed_point import FixedPoint
+from mpi4py import MPI
import args as args
import numpy as np
import time as timer
@@ -21,6 +22,10 @@ tStart = timer.time()
# Parse arguments
arg = args.parse()
+# Who is master process and verbosity
+master = (MPI.COMM_WORLD.Get_rank() == 0)
+verbose = (not arg.quiet) and (master)
+
# Prepare variables of -setnumber for GetDP
setnumber = []
for i in arg.setnumber:
@@ -44,44 +49,46 @@ if(not arg.non_holomorphic):
else:
fixed = FixedPoint(beyn, arg.nhMaxIt, arg.nhMinRratio, not arg.quiet)
# This variable ('fixed') is meant to be called
- if(not arg.quiet): # elsewhere (i.e. Python[] call in GetDP) thanks
+ if(verbose): # elsewhere (i.e. Python[] call in GetDP) thanks
fixed.display() # to the python interpreter.
beyn.display() # Please don't change its name!
l, V, r = fixed.solve()
# Display the computed eigenvalues
-if(not arg.quiet):
+if(verbose):
print
print "Eigenvalues:"
print "----------- "
-for i in range(l.shape[0]):
- if(not arg.quiet): print " # " + str(i) + ": ",
- print "(" + format(np.real(l[i]).tolist(), '+.12e') + ")",
- print "+",
- print "(" + format(np.imag(l[i]).tolist(), '+.12e') + ")*j",
- print "| " + format(r[i], 'e'),
- print "| " + format(r[i]/np.abs(l[i]), 'e')
-
-if(not arg.quiet):
+#if(master):
+# for i in range(l.shape[0]):
+# if(verbose): print " # " + str(i) + ": ",
+# print "(" + format(np.real(l[i]).tolist(), '+.12e') + ")",
+# print "+",
+# print "(" + format(np.imag(l[i]).tolist(), '+.12e') + ")*j",
+# print "| " + format(r[i], 'e'),
+# print "| " + format(r[i]/np.abs(l[i]), 'e')
+#
+if(verbose):
print "----------- "
print "Found: " + str(l.shape[0])
# Generate GetDP views for eigenvectors
-if(not arg.quiet):
+if(verbose):
print
print "Generating eigenvector views..."
-for i in range(l.shape[0]):
- operator.view(V[:, i], "cim" + str(i) + ".pos")
- if(not arg.quiet): print " # View: " + str(i) + " done!"
-
+#if(master):
+# for i in range(l.shape[0]):
+# operator.view(V[:, i], "cim" + str(i) + ".pos")
+# if(not arg.quiet): print " # View: " + str(i) + " done!"
+#
# Stop timer
tStop = timer.time()
# Done
-if(not arg.quiet):
+if(verbose):
print
print "Elapsed time: " + str(round(tStop - tStart)) + "s"
print "Bye!"
diff --git a/bin/fixed_point.py b/bin/fixed_point.py
index d6878a39b73fcfa8bb57e1d36cce16fffab04485..471e40af856489ff4aa9dd13729bf20117e5e786 100644
--- a/bin/fixed_point.py
+++ b/bin/fixed_point.py
@@ -11,8 +11,9 @@ See the LICENSE.txt and README.md for more license and copyright information.
TODO: cluster identification / can a spline be holomorphic?
"""
-import sys as sys
-import numpy as np
+from mpi4py import MPI
+import sys as sys
+import numpy as np
class FixedPoint:
"""A fixed-point solver for handling non-holomorphic operator
@@ -56,7 +57,8 @@ class FixedPoint:
self.__solver = solver
self.__maxIt = maxIt
self.__minNHHRatio = minNHHRatio
- self.__verbose = verbose
+ self.__myProc = MPI.COMM_WORLD.Get_rank()
+ self.__verbose = verbose and (self.__myProc == 0)
# State
self.__cluster = []
diff --git a/bin/solver.py b/bin/solver.py
index 3866d544e96f666f24f9442b70d3cac4fc946dee..834d16802430b934f407ccdf9a19bbdf1d2a6b16 100644
--- a/bin/solver.py
+++ b/bin/solver.py
@@ -7,8 +7,9 @@ Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
-import getdp as GetDP
-import numpy as np
+from mpi4py import MPI # Initializes correctly MPI if solver.py is standalone
+import getdp as GetDP
+import numpy as np
class GetDPWave:
"""A GetDP solver using complex arithmetic for time-harmonic wave problems
diff --git a/example/maxwell_linear/Makefile b/example/maxwell_linear/Makefile
index b17e1e961b382f24e423e11543f5e822692784f5..b711e8189abac8927480a07428308c7fa865ff15 100644
--- a/example/maxwell_linear/Makefile
+++ b/example/maxwell_linear/Makefile
@@ -5,7 +5,7 @@ test:
@echo "make ref: compute the first eigenvalue directly with GetDP"
ref: init
- getdp ref.pro -solve Eig -pos Eig -msh square.msh
+ getdp ref.pro -solve Ref -pos Ref -msh square.msh
cim: init
python ${BIN}/cim.py maxwell.pro square.msh Maxwell 9e8 1e8 -lStart 1 -lStep 1
diff --git a/example/maxwell_linear/maxwell.pro b/example/maxwell_linear/maxwell.pro
index 1e219c873b4763d6cb277ed47d2820411cc47fcf..de19187729d0094195f46ba8aae3268ec7ce1b52 100644
--- a/example/maxwell_linear/maxwell.pro
+++ b/example/maxwell_linear/maxwell.pro
@@ -93,6 +93,7 @@ Resolution{
{ Name A; NameOfFormulation Maxwell; Type Complex; }
}
Operation{
+ MPI_SetCommSelf;
Generate[A];
If(doInit)
diff --git a/example/maxwell_linear/ref.pro b/example/maxwell_linear/ref.pro
index f68e0fa743dec6c475b48f54eb7ced9e3949986b..79518369cf0f186ff56208e94220b4b8c77d9433 100644
--- a/example/maxwell_linear/ref.pro
+++ b/example/maxwell_linear/ref.pro
@@ -56,7 +56,7 @@ FunctionSpace{
}
Formulation{
- { Name Eig; Type FemEquation;
+ { Name Ref; Type FemEquation;
Quantity{
{ Name e; Type Local; NameOfSpace HCurl; }
}
@@ -70,9 +70,9 @@ Formulation{
}
Resolution{
- { Name Eig;
+ { Name Ref;
System{
- { Name A; NameOfFormulation Eig; Type Complex; }
+ { Name A; NameOfFormulation Ref; Type Complex; }
}
Operation{
GenerateSeparate[A];
@@ -82,7 +82,7 @@ Resolution{
}
PostProcessing{
- { Name Eig; NameOfFormulation Eig;
+ { Name Ref; NameOfFormulation Ref;
Quantity{
{ Name e; Value{ Local{ [{e}]; In Omega; Jacobian JVol; } } }
}
@@ -90,7 +90,7 @@ PostProcessing{
}
PostOperation{
- { Name Eig; NameOfPostProcessing Eig;
+ { Name Ref; NameOfPostProcessing Ref;
Operation{
Print[e, OnElementsOf Omega, File "eig.pos", EigenvalueLegend];
}
diff --git a/example/maxwell_sibc/cimResolution.pro b/example/maxwell_sibc/cimResolution.pro
index 99551870563d2775702a773028f245864013abaa..ab7b423c69cd8d052fc2c4448b90c90963bbf388 100644
--- a/example/maxwell_sibc/cimResolution.pro
+++ b/example/maxwell_sibc/cimResolution.pro
@@ -7,6 +7,7 @@ Resolution{
{ Name Solve;
System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
Operation{
+ MPI_SetCommSelf;
Generate[A];
If(doInit)
diff --git a/example/maxwell_sibc_lagrange/cimResolution.pro b/example/maxwell_sibc_lagrange/cimResolution.pro
index 033844d44d1dd5796d1d810e84c61a4b0c116706..186472cc6b6257ff9ebc6a4807a18802ad25bcaa 100644
--- a/example/maxwell_sibc_lagrange/cimResolution.pro
+++ b/example/maxwell_sibc_lagrange/cimResolution.pro
@@ -7,6 +7,7 @@ Resolution{
{ Name Solve;
System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
Operation{
+ MPI_SetCommSelf;
If(!doInit)
Evaluate[Python[]{"real_corrector_init.py"}];
EndIf