diff --git a/Examples/demo_2D_real-time.py b/Examples/demo_2D_real-time.py
index 6b6f28df292644995174cf41e35d87121854494e..4ed123c8035f71b87ebc9ef07f9e5bdddb162a45 100755
--- a/Examples/demo_2D_real-time.py
+++ b/Examples/demo_2D_real-time.py
@@ -1,29 +1,28 @@
#!/usr/bin/env python
-import time
import parmepy
from parmepy.domain.box import Box
from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU
+#parmepy.gpu.CL_PROFILE = True
+
from parmepy.fields.continuous import Field
from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem
from parmepy.operator.analytic import Analytic
from parmepy.gpu.QtRendering import QtOpenGLRendering
+from parmepy.problem.simulation import Simulation
import math
-import sys
import numpy as np
-parmepy.gpu.CL_PROFILE = False
-
norm2 = lambda x, y: x * x + y * y
norm1 = lambda x, y: abs(x) + abs(y)
norminf = lambda x, y: max(abs(x), abs(y))
def initScalar(x, y):
- return math.exp(-(norm2(x - 0.5, y - 0.75)/0.0225)**6) \
- + 0.75*math.exp(-(norm2(x - 0.75, y - 0.25)/0.0225)**6) \
- + 0.5*math.exp(-(norm1(x - 0.4, y - 0.4)/0.1)**6) \
- + 0.25*math.exp(-(norminf(x - 0.6, y - 0.5)/0.08)**6)
+ return math.exp(-(norm2(x - 0.5, y - 0.75) / 0.0225) ** 6) \
+ + 0.75 * math.exp(-(norm2(x - 0.75, y - 0.25) / 0.0225) ** 6) \
+ + 0.5 * math.exp(-(norm1(x - 0.4, y - 0.4) / 0.1) ** 6) \
+ + 0.25 * math.exp(-(norminf(x - 0.6, y - 0.5) / 0.08) ** 6)
# if norm2(x - 0.5, y - 0.75) < 0.0225:
# return 0.25
# elif norm2(x - 0.75, y - 0.25) < 0.0225:
@@ -36,62 +35,66 @@ def initScalar(x, y):
# return 0.
-def run(nb=257):
- dim = 2
- boxLength = [1., 1.]
- boxMin = [0., 0.]
- if isinstance(nb, list):
- nbElem = nb
- else:
- nbElem = [nb, nb]
-
- timeStep = 0.075
- finalTime = 3.0 + timeStep
- ## Domain
- box = Box(dim, length=boxLength, origin=boxMin)
-
- ## Fields
- scal = Field(domain=box, name='Scalar', formula=initScalar)
- velo = Field(domain=box, name='Velocity', isVector=True)
-
- ## Operators
- advec = Advection(velo, scal,
- resolutions={velo: nbElem,
- scal: nbElem},
- method='gpu_1k_m6prime',
- #src=['./demo_2D.cl'],
- precision=PARMES_REAL_GPU,
- splittingConfig='o2'
- )
- velocity = Analytic(velo,
- resolutions={velo: nbElem},
- method='gpu',
- src=['./demo_2D.cl'],
- precision=PARMES_REAL_GPU,
- )
- render = QtOpenGLRendering(scal)
-
- # Problem
- pb = TransportProblem([velocity, advec],
- monitors=[render])
-
- # Setting solver to Problem
- pb.setUp(finalTime, timeStep)
-
- # We copy the first discretisation of scalar
- scalar_initial = np.copy(scal.discreteFields.values()[0].data)
-
- ## Solve problem
- pb.solve()
-
- scal_disc = scal.discreteFields.values()[0]
- scal_disc.toHost()
- print 'Erreur : ', np.max(np.abs(scalar_initial -
- scal_disc.data)) / np.max(scalar_initial)
- print np.linalg.norm(scalar_initial - scal_disc.data, ord=2)/np.linalg.norm(scalar_initial, ord=2)
-
- pb.finalize()
-
-if __name__ == "__main__":
- run([513, 513])
-
+def vitesse(x, y, t=0):
+ vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2.) * \
+ math.cos(t * math.pi / 3.)
+ vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2.) * \
+ math.cos(t * math.pi / 3.)
+ return vx, vy
+
+
+dim = 2
+boxLength = [1., 1.]
+boxMin = [0., 0.]
+nbElem = [513, 513]
+
+timeStep = 0.075
+finalTime = 3.0 + timeStep
+
+simu = Simulation(0., finalTime, int(finalTime / timeStep))
+
+## Domain
+box = Box(dim, length=boxLength, origin=boxMin)
+
+## Fields
+scal = Field(domain=box, name='Scalar', formula=initScalar)
+velo = Field(domain=box, name='Velocity', isVector=True)
+
+## Operators
+advec = Advection(velo, scal,
+ resolutions={velo: nbElem,
+ scal: nbElem},
+ method='gpu_1k_m6prime',
+ #src=['./demo_2D.cl'],
+ precision=PARMES_REAL_GPU,
+ splittingConfig='o2'
+ )
+velocity = Analytic(velo, # formula=vitesse,
+ resolutions={velo: nbElem},
+ method='gpu',
+ src=['./demo_2D.cl'],
+ precision=PARMES_REAL_GPU,
+ )
+render = QtOpenGLRendering(scal)
+
+# Problem
+pb = TransportProblem([velocity, advec], simu,
+ monitors=[render])
+
+# Setting solver to Problem
+pb.setUp()
+
+# We copy the first discretisation of scalar
+scalar_initial = np.copy(scal.discreteFields.values()[0].data)
+
+## Solve problem
+pb.solve()
+
+scal_disc = scal.discreteFields.values()[0]
+scal_disc.toHost()
+print 'Erreur : ', np.max(np.abs(scalar_initial -
+ scal_disc.data)) / np.max(scalar_initial)
+print np.linalg.norm(scalar_initial - scal_disc.data, ord=2) / \
+ np.linalg.norm(scalar_initial, ord=2)
+
+pb.finalize()
diff --git a/Examples/levelSet2D.py b/Examples/levelSet2D.py
index a0a2e0de9e71f9b33117a1c032a5c22cc9e54b68..84a59b5f96a42b5a020caaeb3c7916be6994e290 100755
--- a/Examples/levelSet2D.py
+++ b/Examples/levelSet2D.py
@@ -1,119 +1,99 @@
#!/usr/bin/env python
-import time
from parmepy.domain.box import Box
from parmepy.gpu import PARMES_REAL_GPU
from parmepy.fields.continuous import Field
from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem
from parmepy.operator.monitors.printer import Printer
+from parmepy.problem.simulation import Simulation
from parmepy.operator.analytic import Analytic
-import math
-import sys
-
+# import math
# def vitesse(x, y):
# vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2)
# vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2)
# return vx, vy
-def vitesse(x, y, t, dt, ite):
- vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2)*math.cos(t*math.pi/3.0)
- vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2)*math.cos(t*math.pi/3.0)
- return vx, vy
-
-
-def scalaire(x, y):
- rr = math.sqrt((x - 0.5) ** 2 + (y - 0.75) ** 2)
- if rr < 0.15:
- return 1.
- else:
- return 0.
-
-
-def run(nb=257):
- dim = 2
- boxLength = [1., 1.]
- boxMin = [0., 0.]
- if isinstance(nb, list):
- nbElem = nb
- else:
- nbElem = [nb, nb]
-
- timeStep = 0.025
- finalTime = 3.
- outputFilePrefix = './res_2D/levelSet_2D_rect_'
- outputModulo = 1
-
- t0 = time.time()
-
- ## Domain
- box = Box(dim, length=boxLength, origin=boxMin)
-
- ## Fields
- scal = Field(domain=box, name='Scalar')
- velo = Field(domain=box, name='Velocity', isVector=True)
-
- ## Operators
- advec = Advection(velo, scal,
- resolutions={velo: nbElem,
- scal: nbElem},
- method='gpu_1k_m4prime',
- #method='gpu_1k_m6prime',
- #method='gpu_1k_m8prime',
- #method='gpu_2k_m4prime',
- #method='gpu_2k_m6prime',
- #method='gpu_2k_m8prime',
- #method='scales'
- src=['./levelSet2D.cl'],
- precision=PARMES_REAL_GPU,
- #precision=PARMES_DOUBLE_GPU,
- #splittingConfig='o2'
- #splittingConfig='y_only'
- splittingConfig='o2_FullHalf'
- )
- velocity = Analytic(velo,
- resolutions={velo: nbElem},
- method='gpu',
- src=['./levelSet2D.cl'],
- precision=PARMES_REAL_GPU,
- )
-
- ##Problem
- #pb = TransportProblem([advec],
- pb = TransportProblem([velocity, advec],
- monitors=[Printer(fields=[scal],
- frequency=outputModulo,
- outputPrefix=outputFilePrefix)])
-
- ## Setting solver to Problem
- pb.setUp(finalTime, timeStep)
-
- t1 = time.time()
- ## Solve problem
- timings = pb.solve()
- tf = time.time()
-
- print "\n"
- print "Total time : ", tf - t0, "sec (CPU)"
- print "Init time : ", t1 - t0, "sec (CPU)"
- print "Solving time : ", tf - t1, "sec (CPU)"
- print ""
- return (pb.time_info[0] + "\"Solving time\" \"Init time\" \"Total time\"",
- pb.time_info[1] + str(tf - t1) + " " + str(t1 - t0) + " " + str(tf - t0))
-
-if __name__ == "__main__":
- run(257)
- # timings = {}
- # f = open('bench_levelSet2D.dat', 'w')
- # header = "#size dim nPart "
- # sizes = xrange(128, 4096+2048+1, 128)
- # sizes = xrange(128, 512+1, 128)
- # sizes = [1024,2048,4096,4096+2048]
- # for s in sizes:
- # h, t = run(s)
- # timings[s] = t
- # header += h
- # f.write(header + "\n")
- # for s in sizes:
- # f.write(str(s) + " " + str(2) + " " + str(s**2) + " " + str(timings[s]) + "\n")
- # f.close()
+
+# def vitesse(x, y, t=0):
+# vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2) * \
+# math.cos(t * math.pi / 3.)
+# vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2) * \
+# math.cos(t * math.pi / 3.)
+# return vx, vy
+
+
+# def scalaire(x, y):
+# rr = math.sqrt((x - 0.5) ** 2 + (y - 0.75) ** 2)
+# if rr < 0.15:
+# return 1.
+# else:
+# return 0.
+
+
+dim = 2
+boxLength = [1., 1.]
+boxMin = [0., 0.]
+nbElem = [1025, 1025]
+
+timeStep = 0.025
+finalTime = 3.
+outputFilePrefix = './res_2D/levelSet_2D_rect_'
+outputModulo = 0
+simu = Simulation(0., finalTime, int(finalTime / timeStep))
+
+## Domain
+box = Box(dim, length=boxLength, origin=boxMin)
+
+## Fields
+scal = Field(domain=box, name='Scalar')
+velo = Field(domain=box, name='Velocity', isVector=True)
+
+## Operators
+advec = Advection(velo, scal,
+ resolutions={velo: nbElem,
+ scal: nbElem},
+ method='gpu_2k_m4prime',
+ #method='gpu_1k_m6prime',
+ #method='gpu_1k_m8prime',
+ #method='gpu_2k_m4prime',
+ #method='gpu_2k_m6prime',
+ #method='gpu_2k_m8prime',
+ #method='scales'
+ src=['./levelSet2D.cl'],
+ precision=PARMES_REAL_GPU,
+ #precision=PARMES_DOUBLE_GPU,
+ #splittingConfig='o2'
+ #splittingConfig='y_only'
+ #splittingConfig='o2_FullHalf'
+ )
+velocity = Analytic(velo,
+ resolutions={velo: nbElem},
+ method='gpu',
+ src=['./levelSet2D.cl'],
+ precision=PARMES_REAL_GPU,
+ )
+
+##Problem
+#pb = TransportProblem([advec],simu,
+pb = TransportProblem([velocity, advec], simu,
+ monitors=[Printer(fields=[scal],
+ frequency=outputModulo,
+ prefix=outputFilePrefix)])
+
+## Setting solver to Problem
+pb.setUp()
+
+## Solve problem
+pb.solve()
+
+pb.finalize()
+
+th, tt = pb.timer.toString()
+for op in pb.operators:
+ oth, ott = op.timer.toString()
+ th += oth
+ tt += ott
+print th
+print tt
+
diff --git a/Examples/levelSet3D.py b/Examples/levelSet3D.py
index 9a78e45696a666b96c9398302755fe41930b5367..23afc96eb98d1b7cb6e7d3064b1f12275087e00f 100755
--- a/Examples/levelSet3D.py
+++ b/Examples/levelSet3D.py
@@ -1,5 +1,4 @@
#!/usr/bin/env python
-import time
import parmepy as pp
from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU
from parmepy.fields.continuous import Field
@@ -7,100 +6,81 @@ from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem
from parmepy.operator.monitors.printer import Printer
from parmepy.operator.analytic import Analytic
-from math import sin, cos, pi
-import sys
-
-
-def vitesse(x, y, z, t=0, dt=0, ite=0):
- vx = 2. * sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z)*cos(t*pi/3.)
- vy = -sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)*cos(t*pi/3.)
- vz = -sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2*cos(t*pi/3.)
- return vx, vy, vz
-
-
-def run(nb=257):
- dim = 3
- boxLength = [1., 1., 1.]
- boxMin = [0., 0., 0.]
- if isinstance(nb, list):
- nbElem = nb
- else:
- nbElem = [nb, nb, nb]
-
- timeStep = 0.05
- finalTime = 3.
- outputFilePrefix = './res3D_new/levelSet_3D_'
- outputModulo = 1
-
- t0 = time.time()
-
- ## Domain
- box = pp.Box(dim, length=boxLength, origin=boxMin)
-
- ## Fields
- scal = Field(domain=box, name='Scalar')
- velo = Field(domain=box, name='Velocity', isVector=True)
-
- ## Operators
- advec = Advection(velo, scal,
- resolutions={velo: nbElem,
- scal: nbElem},
- #method='scales',
- method='gpu_1k_m4prime',
- #method='gpu_1k_m6prime',
- #method='gpu_1k_m8prime',
- #method='gpu_2k_m4prime',
- #method='gpu_2k_m6prime',
- #method='gpu_2k_m8prime',
- src=['./levelSet3D.cl'],
- precision=PARMES_REAL_GPU,
- #precision=PARMES_DOUBLE_GPU,
- splittingConfig='o2'
- #splittingConfig='o2_FullHalf'
- )
- velocity = Analytic(velo, # formula=vitesse,
- resolutions={velo: nbElem},
- method='gpu',
- src=['./levelSet3D.cl'],
- precision=PARMES_REAL_GPU,
- )
-
- ##Problem
- #pb = TransportProblem([advec],
- pb = TransportProblem([velocity, advec],
- monitors=[Printer(fields=[scal],
- frequency=outputModulo,
- outputPrefix=outputFilePrefix)])
-
- ## Setting solver to Problem
- pb.setUp(finalTime, timeStep)
-
- t1 = time.time()
-
- ## Solve problem
- timings = pb.solve()
-
- tf = time.time()
-
- print "\n"
- print "Total time : ", tf - t0, "sec (CPU)"
- print "Init time : ", t1 - t0, "sec (CPU)"
- print "Solving time : ", tf - t1, "sec (CPU)"
- print ""
- return pb.time_info[0] + "\"Solving time\" \"Init time\" \"Total time\"", pb.time_info[1] + str(tf - t1) + " " + str(t1 - t0) + " " + str(tf - t0)
-
-if __name__ == "__main__":
- run([129, 129, 129])
- # timings = {}
- # f = open('bench_levelSet3D.dat', 'w')
- # header = "#size dim nPart "
- # sizes = xrange(32, 256+32*3+1, 32)
- # #sizes = [256+32+32+32]
- # for s in sizes:
- # h, t = run(s)
- # timings[s] = t
- # header += h
- # f.write(header + "\n")
- # for s in sizes:
- # f.write(str(s) + " " + str(3) + " " + str(s**3) + " " + str(timings[s]) + "\n")
- # f.close()
+from parmepy.problem.simulation import Simulation
+# from math import sin, cos, pi
+
+
+# def scalaire(x,y,z):
+# return 1.
+
+
+# def vitesse(x, y, z, t=0):
+# vx = 2. * sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z)*cos(t*pi/3.)
+# vy = -sin(2.*pi*x)*sin(pi*y)**2*sin(2.*pi*z)*cos(t*pi/3.)
+# vz = -sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2*cos(t*pi/3.)
+# return vx, vy, vz
+
+dim = 3
+boxLength = [1., 1., 1.]
+boxMin = [0., 0., 0.]
+nbElem = 3 * [129]
+
+timeStep = 0.05
+finalTime = 3.
+outputFilePrefix = './res3D_new/levelSet_3D_'
+outputModulo = 1
+simu = Simulation(0., finalTime, int(finalTime / timeStep))
+
+## Domain
+box = pp.Box(dim, length=boxLength, origin=boxMin)
+
+## Fields
+scal = Field(domain=box, name='Scalar')
+velo = Field(domain=box, name='Velocity', isVector=True)
+
+## Operators
+advec = Advection(velo, scal,
+ resolutions={velo: nbElem,
+ scal: nbElem},
+ #method='scales',
+ method='gpu_1k_m4prime',
+ #method='gpu_1k_m6prime',
+ #method='gpu_1k_m8prime',
+ #method='gpu_2k_m4prime',
+ #method='gpu_2k_m6prime',
+ #method='gpu_2k_m8prime',
+ src=['./levelSet3D.cl'],
+ precision=PARMES_REAL_GPU,
+ #precision=PARMES_DOUBLE_GPU,
+ splittingConfig='o2'
+ #splittingConfig='o2_FullHalf'
+ )
+velocity = Analytic(velo,
+ resolutions={velo: nbElem},
+ method='gpu',
+ src=['./levelSet3D.cl'],
+ precision=PARMES_REAL_GPU,
+ )
+
+##Problem
+#pb = TransportProblem([advec],
+pb = TransportProblem([velocity, advec], simu,
+ monitors=[Printer(fields=[scal],
+ frequency=outputModulo,
+ prefix=outputFilePrefix)])
+
+## Setting solver to Problem
+pb.setUp()
+
+## Solve problem
+pb.solve()
+
+pb.finalize()
+
+th, tt = pb.timer.toString()
+for op in pb.operators:
+ oth, ott = op.timer.toString()
+ th += oth
+ tt += ott
+print th
+print tt