Skip to content
Snippets Groups Projects
Commit 27a862a9 authored by Jean-Matthieu Etancelin's avatar Jean-Matthieu Etancelin
Browse files

Update examples

parent 7049a2d5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Examples/demo_2D_real-time.py
+ 70
− 67
View file @ 27a862a9
#!/usr/bin/env python #!/usr/bin/env python
import time
import parmepy import parmepy
from parmepy.domain.box import Box from parmepy.domain.box import Box
from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU
#parmepy.gpu.CL_PROFILE = True
from parmepy.fields.continuous import Field from parmepy.fields.continuous import Field
from parmepy.operator.advection import Advection from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem from parmepy.problem.transport import TransportProblem
from parmepy.operator.analytic import Analytic from parmepy.operator.analytic import Analytic
from parmepy.gpu.QtRendering import QtOpenGLRendering from parmepy.gpu.QtRendering import QtOpenGLRendering
from parmepy.problem.simulation import Simulation
import math import math
import sys
import numpy as np import numpy as np
parmepy.gpu.CL_PROFILE = False
norm2 = lambda x, y: x * x + y * y norm2 = lambda x, y: x * x + y * y
norm1 = lambda x, y: abs(x) + abs(y) norm1 = lambda x, y: abs(x) + abs(y)
norminf = lambda x, y: max(abs(x), abs(y)) norminf = lambda x, y: max(abs(x), abs(y))
def initScalar(x, y): def initScalar(x, y):
return math.exp(-(norm2(x - 0.5, y - 0.75)/0.0225)**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.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.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) + 0.25 * math.exp(-(norminf(x - 0.6, y - 0.5) / 0.08) ** 6)
# if norm2(x - 0.5, y - 0.75) < 0.0225: # if norm2(x - 0.5, y - 0.75) < 0.0225:
# return 0.25 # return 0.25
# elif norm2(x - 0.75, y - 0.25) < 0.0225: # elif norm2(x - 0.75, y - 0.25) < 0.0225:
...@@ -36,62 +35,66 @@ def initScalar(x, y): ...@@ -36,62 +35,66 @@ def initScalar(x, y):
# return 0. # return 0.
def run(nb=257): def vitesse(x, y, t=0):
dim = 2 vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2.) * \
boxLength = [1., 1.] math.cos(t * math.pi / 3.)
boxMin = [0., 0.] vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2.) * \
if isinstance(nb, list): math.cos(t * math.pi / 3.)
nbElem = nb return vx, vy
else:
nbElem = [nb, nb]
dim = 2
timeStep = 0.075 boxLength = [1., 1.]
finalTime = 3.0 + timeStep boxMin = [0., 0.]
## Domain nbElem = [513, 513]
box = Box(dim, length=boxLength, origin=boxMin)
timeStep = 0.075
## Fields finalTime = 3.0 + timeStep
scal = Field(domain=box, name='Scalar', formula=initScalar)
velo = Field(domain=box, name='Velocity', isVector=True) simu = Simulation(0., finalTime, int(finalTime / timeStep))
## Operators ## Domain
advec = Advection(velo, scal, box = Box(dim, length=boxLength, origin=boxMin)
resolutions={velo: nbElem,
scal: nbElem}, ## Fields
method='gpu_1k_m6prime', scal = Field(domain=box, name='Scalar', formula=initScalar)
#src=['./demo_2D.cl'], velo = Field(domain=box, name='Velocity', isVector=True)
precision=PARMES_REAL_GPU,
splittingConfig='o2' ## Operators
) advec = Advection(velo, scal,
velocity = Analytic(velo, resolutions={velo: nbElem,
resolutions={velo: nbElem}, scal: nbElem},
method='gpu', method='gpu_1k_m6prime',
src=['./demo_2D.cl'], #src=['./demo_2D.cl'],
precision=PARMES_REAL_GPU, precision=PARMES_REAL_GPU,
) splittingConfig='o2'
render = QtOpenGLRendering(scal) )
velocity = Analytic(velo, # formula=vitesse,
# Problem resolutions={velo: nbElem},
pb = TransportProblem([velocity, advec], method='gpu',
monitors=[render]) src=['./demo_2D.cl'],
precision=PARMES_REAL_GPU,
# Setting solver to Problem )
pb.setUp(finalTime, timeStep) render = QtOpenGLRendering(scal)
# We copy the first discretisation of scalar # Problem
scalar_initial = np.copy(scal.discreteFields.values()[0].data) pb = TransportProblem([velocity, advec], simu,
monitors=[render])
## Solve problem
pb.solve() # Setting solver to Problem
pb.setUp()
scal_disc = scal.discreteFields.values()[0]
scal_disc.toHost() # We copy the first discretisation of scalar
print 'Erreur : ', np.max(np.abs(scalar_initial - scalar_initial = np.copy(scal.discreteFields.values()[0].data)
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) ## Solve problem
pb.solve()
pb.finalize()
scal_disc = scal.discreteFields.values()[0]
if __name__ == "__main__": scal_disc.toHost()
run([513, 513]) 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()
Examples/levelSet2D.py
+ 85
− 105
View file @ 27a862a9
#!/usr/bin/env python #!/usr/bin/env python
import time
from parmepy.domain.box import Box from parmepy.domain.box import Box
from parmepy.gpu import PARMES_REAL_GPU from parmepy.gpu import PARMES_REAL_GPU
from parmepy.fields.continuous import Field from parmepy.fields.continuous import Field
from parmepy.operator.advection import Advection from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem from parmepy.problem.transport import TransportProblem
from parmepy.operator.monitors.printer import Printer from parmepy.operator.monitors.printer import Printer
from parmepy.problem.simulation import Simulation
from parmepy.operator.analytic import Analytic from parmepy.operator.analytic import Analytic
import math # import math
import sys
# def vitesse(x, y): # def vitesse(x, y):
# vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2) # 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) # vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2)
# return vx, vy # 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) # def vitesse(x, y, t=0):
vy = math.sin(y * math.pi) ** 2 * math.sin(x * math.pi * 2)*math.cos(t*math.pi/3.0) # vx = -math.sin(x * math.pi) ** 2 * math.sin(y * math.pi * 2) * \
return vx, vy # 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.)
def scalaire(x, y): # return vx, vy
rr = math.sqrt((x - 0.5) ** 2 + (y - 0.75) ** 2)
if rr < 0.15:
return 1. # def scalaire(x, y):
else: # rr = math.sqrt((x - 0.5) ** 2 + (y - 0.75) ** 2)
return 0. # if rr < 0.15:
# return 1.
# else:
def run(nb=257): # return 0.
dim = 2
boxLength = [1., 1.]
boxMin = [0., 0.] dim = 2
if isinstance(nb, list): boxLength = [1., 1.]
nbElem = nb boxMin = [0., 0.]
else: nbElem = [1025, 1025]
nbElem = [nb, nb]
timeStep = 0.025
timeStep = 0.025 finalTime = 3.
finalTime = 3. outputFilePrefix = './res_2D/levelSet_2D_rect_'
outputFilePrefix = './res_2D/levelSet_2D_rect_' outputModulo = 0
outputModulo = 1 simu = Simulation(0., finalTime, int(finalTime / timeStep))
t0 = time.time() ## Domain
box = Box(dim, length=boxLength, origin=boxMin)
## Domain
box = Box(dim, length=boxLength, origin=boxMin) ## Fields
scal = Field(domain=box, name='Scalar')
## Fields velo = Field(domain=box, name='Velocity', isVector=True)
scal = Field(domain=box, name='Scalar')
velo = Field(domain=box, name='Velocity', isVector=True) ## Operators
advec = Advection(velo, scal,
## Operators resolutions={velo: nbElem,
advec = Advection(velo, scal, scal: nbElem},
resolutions={velo: nbElem, method='gpu_2k_m4prime',
scal: nbElem}, #method='gpu_1k_m6prime',
method='gpu_1k_m4prime', #method='gpu_1k_m8prime',
#method='gpu_1k_m6prime', #method='gpu_2k_m4prime',
#method='gpu_1k_m8prime', #method='gpu_2k_m6prime',
#method='gpu_2k_m4prime', #method='gpu_2k_m8prime',
#method='gpu_2k_m6prime', #method='scales'
#method='gpu_2k_m8prime', src=['./levelSet2D.cl'],
#method='scales' precision=PARMES_REAL_GPU,
src=['./levelSet2D.cl'], #precision=PARMES_DOUBLE_GPU,
precision=PARMES_REAL_GPU, #splittingConfig='o2'
#precision=PARMES_DOUBLE_GPU, #splittingConfig='y_only'
#splittingConfig='o2' #splittingConfig='o2_FullHalf'
#splittingConfig='y_only' )
splittingConfig='o2_FullHalf' velocity = Analytic(velo,
) resolutions={velo: nbElem},
velocity = Analytic(velo, method='gpu',
resolutions={velo: nbElem}, src=['./levelSet2D.cl'],
method='gpu', precision=PARMES_REAL_GPU,
src=['./levelSet2D.cl'], )
precision=PARMES_REAL_GPU,
) ##Problem
#pb = TransportProblem([advec],simu,
##Problem pb = TransportProblem([velocity, advec], simu,
#pb = TransportProblem([advec], monitors=[Printer(fields=[scal],
pb = TransportProblem([velocity, advec], frequency=outputModulo,
monitors=[Printer(fields=[scal], prefix=outputFilePrefix)])
frequency=outputModulo,
outputPrefix=outputFilePrefix)]) ## Setting solver to Problem
pb.setUp()
## Setting solver to Problem
pb.setUp(finalTime, timeStep) ## Solve problem
pb.solve()
t1 = time.time()
## Solve problem pb.finalize()
timings = pb.solve()
tf = time.time() th, tt = pb.timer.toString()
for op in pb.operators:
print "\n" oth, ott = op.timer.toString()
print "Total time : ", tf - t0, "sec (CPU)" th += oth
print "Init time : ", t1 - t0, "sec (CPU)" tt += ott
print "Solving time : ", tf - t1, "sec (CPU)" print th
print "" print tt
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()
Examples/levelSet3D.py
+ 78
− 98
View file @ 27a862a9
#!/usr/bin/env python #!/usr/bin/env python
import time
import parmepy as pp import parmepy as pp
from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU from parmepy.gpu import PARMES_REAL_GPU, PARMES_DOUBLE_GPU
from parmepy.fields.continuous import Field from parmepy.fields.continuous import Field
...@@ -7,100 +6,81 @@ from parmepy.operator.advection import Advection ...@@ -7,100 +6,81 @@ from parmepy.operator.advection import Advection
from parmepy.problem.transport import TransportProblem from parmepy.problem.transport import TransportProblem
from parmepy.operator.monitors.printer import Printer from parmepy.operator.monitors.printer import Printer
from parmepy.operator.analytic import Analytic from parmepy.operator.analytic import Analytic
from math import sin, cos, pi from parmepy.problem.simulation import Simulation
import sys # from math import sin, cos, pi
def vitesse(x, y, z, t=0, dt=0, ite=0): # def scalaire(x,y,z):
vx = 2. * sin(pi*x)**2*sin(2.*pi*y)*sin(2.*pi*z)*cos(t*pi/3.) # return 1.
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 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.)
def run(nb=257): # vz = -sin(2.*pi*x)*sin(2.*pi*y)*sin(pi*z)**2*cos(t*pi/3.)
dim = 3 # return vx, vy, vz
boxLength = [1., 1., 1.]
boxMin = [0., 0., 0.] dim = 3
if isinstance(nb, list): boxLength = [1., 1., 1.]
nbElem = nb boxMin = [0., 0., 0.]
else: nbElem = 3 * [129]
nbElem = [nb, nb, nb]
timeStep = 0.05
timeStep = 0.05 finalTime = 3.
finalTime = 3. outputFilePrefix = './res3D_new/levelSet_3D_'
outputFilePrefix = './res3D_new/levelSet_3D_' outputModulo = 1
outputModulo = 1 simu = Simulation(0., finalTime, int(finalTime / timeStep))
t0 = time.time() ## Domain
box = pp.Box(dim, length=boxLength, origin=boxMin)
## Domain
box = pp.Box(dim, length=boxLength, origin=boxMin) ## Fields
scal = Field(domain=box, name='Scalar')
## Fields velo = Field(domain=box, name='Velocity', isVector=True)
scal = Field(domain=box, name='Scalar')
velo = Field(domain=box, name='Velocity', isVector=True) ## Operators
advec = Advection(velo, scal,
## Operators resolutions={velo: nbElem,
advec = Advection(velo, scal, scal: nbElem},
resolutions={velo: nbElem, #method='scales',
scal: nbElem}, method='gpu_1k_m4prime',
#method='scales', #method='gpu_1k_m6prime',
method='gpu_1k_m4prime', #method='gpu_1k_m8prime',
#method='gpu_1k_m6prime', #method='gpu_2k_m4prime',
#method='gpu_1k_m8prime', #method='gpu_2k_m6prime',
#method='gpu_2k_m4prime', #method='gpu_2k_m8prime',
#method='gpu_2k_m6prime', src=['./levelSet3D.cl'],
#method='gpu_2k_m8prime', precision=PARMES_REAL_GPU,
src=['./levelSet3D.cl'], #precision=PARMES_DOUBLE_GPU,
precision=PARMES_REAL_GPU, splittingConfig='o2'
#precision=PARMES_DOUBLE_GPU, #splittingConfig='o2_FullHalf'
splittingConfig='o2' )
#splittingConfig='o2_FullHalf' velocity = Analytic(velo,
) resolutions={velo: nbElem},
velocity = Analytic(velo, # formula=vitesse, method='gpu',
resolutions={velo: nbElem}, src=['./levelSet3D.cl'],
method='gpu', precision=PARMES_REAL_GPU,
src=['./levelSet3D.cl'], )
precision=PARMES_REAL_GPU,
) ##Problem
#pb = TransportProblem([advec],
##Problem pb = TransportProblem([velocity, advec], simu,
#pb = TransportProblem([advec], monitors=[Printer(fields=[scal],
pb = TransportProblem([velocity, advec], frequency=outputModulo,
monitors=[Printer(fields=[scal], prefix=outputFilePrefix)])
frequency=outputModulo,
outputPrefix=outputFilePrefix)]) ## Setting solver to Problem
pb.setUp()
## Setting solver to Problem
pb.setUp(finalTime, timeStep) ## Solve problem
pb.solve()
t1 = time.time()
pb.finalize()
## Solve problem
timings = pb.solve() th, tt = pb.timer.toString()
for op in pb.operators:
tf = time.time() oth, ott = op.timer.toString()
th += oth
print "\n" tt += ott
print "Total time : ", tf - t0, "sec (CPU)" print th
print "Init time : ", t1 - t0, "sec (CPU)" print tt
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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment