suzanne.soy/FreeCAD
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Developed time integrator for fluid

Browse Source
This commit is contained in:
Jose Luis Cercós Pita 2012-08-04 17:11:41 +02:00 committed by Yorik van Havre
parent 9bc0f7eea4
commit bc181f3df3
6 changed files with 154 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/Mod/Ship/CMakeLists.txt
View File

@ -152,6 +152,7 @@ SET(SimRun_SRCS
simRun/Sim/initialization.py
simRun/Sim/matrixGen.py
simRun/Sim/computeSources.py
simRun/Sim/fsEvolution.py
)
SOURCE_GROUP("simrun" FILES ${SimRun_SRCS})

3
src/Mod/Ship/Makefile.am
View File

@ -104,7 +104,8 @@ nobase_data_DATA = \
simRun/Sim/__init__.py \
simRun/Sim/initialization.py \
simRun/Sim/matrixGen.py \
simRun/Sim/computeSources.py
simRun/Sim/computeSources.py \
simRun/Sim/fsEvolution.py
CLEANFILES = $(BUILT_SOURCES)

1
src/Mod/Ship/simRun/Sim/__init__.py
View File

@ -24,3 +24,4 @@
from initialization import *
from matrixGen import *
from computeSources import *
from fsEvolution import *

136
src/Mod/Ship/simRun/Sim/fsEvolution.py Normal file
View File

@ -0,0 +1,136 @@
#***************************************************************************
#* *
#* Copyright (c) 2011, 2012 *
#* Jose Luis Cercos Pita <jlcercos@gmail.com> *
#* *
#* This program is free software; you can redistribute it and/or modify *
#* it under the terms of the GNU Lesser General Public License (LGPL) *
#* as published by the Free Software Foundation; either version 2 of *
#* the License, or (at your option) any later version. *
#* for detail see the LICENCE text file. *
#* *
#* 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 Library General Public License for more details. *
#* *
#* You should have received a copy of the GNU Library General Public *
#* License along with this program; if not, write to the Free Software *
#* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *
#* USA *
#* *
#***************************************************************************
# numpy
import numpy as np
grav=9.81
class simFSEvolution:
def __init__(self, context=None, queue=None):
""" Constructor.
@param context OpenCL context where apply. Only for compatibility,
must be None.
@param queue OpenCL command queue. Only for compatibility,
must be None.
"""
self.context = context
self.queue = queue
def execute(self, fs, waves, dt, t):
""" Compute free surface for next time step.
@param fs Free surface instance.
@param waves Waves instance.
@param dt Time step.
@param t Actual time (without adding dt).
"""
self.fs = fs
# Allocate memory
nx = self.fs['Nx']
ny = self.fs['Ny']
nF = nx*ny
# Evaluate potential gradients
grad = self.evaluateGradient()
# Integrate variables
for i in range(0,nx):
for j in range(0,ny):
# Free surface points position
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*grad[i*ny+j][2]
# Velocity potential
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \
dt*self.fs['accPot'][i,j] - \
0.5*dt*dt*grav*grad[i*ny+j][2]
# Acceleration potential
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - \
dt*grav*grad[i*ny+j][2]
# Force boundary conditions
for i in range(0,nx):
for j in [0,ny-1]:
self.boundaryCondition(i,j, waves, dt, t)
def evaluateGradient(self):
""" Evaluate potential gradients over free surface.
@return Potential gradients.
"""
nx = self.fs['Nx']
ny = self.fs['Ny']
nF = nx*ny
grad = np.ndarray((nF,3), dtype=np.float32)
for i in range(0,nx):
for j in range(0,ny):
pos = self.fs['pos'][i,j]
grad[i*ny+j] = self.gradientphi(pos)
return grad
def gradientphi(self, pos):
""" Compute gradient over desired position.
@param pos Point to evaluate.
@return Potential gradient.
"""
nx = self.fs['Nx']
ny = self.fs['Ny']
grad = np.ndarray(3, dtype=np.float32)
for i in range(0,nx):
for j in range(0,ny):
# Get source position (desingularized)
srcPos = np.copy(self.fs['pos'][i,j])
area = self.fs['area'][i,j]
srcPos[2] = srcPos[2] + np.sqrt(area)
src = self.fs['velSrc'][i,j]
# Get distance between points
d = pos-srcPos
grad = grad + d/np.dot(d,d)*src*area
return grad
def boundaryCondition(self, i,j, waves, dt, t):
""" Compute free surface at boundaries, assuming that only
incident wave can be taken into account.
@param i First free surface cell index.
@param j Second free surface cell index.
@param waves Waves instance.
@param dt Time step.
@param t Actual time (without adding dt).
"""
pos = self.fs['pos'][i,j]
pos[2] = 0.
for w in waves['data']:
A = w[0]
T = w[1]
phase = w[2]
heading = np.pi*w[3]/180.0
wl = 0.5 * grav / np.pi * T*T
k = 2.0*np.pi/wl
frec = 2.0*np.pi/T
pos = self.fs['pos'][i,j]
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
amp = A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + amp
amp = frec*A*np.cos(k*l - frec*(t+dt) + phase)
self.fs['vel'][i,j][2] = self.fs['vel'][i,j][2] - amp
amp = frec*frec*A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['acc'][i,j][2] = self.fs['acc'][i,j][2] - amp
amp = grav/frec*A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + amp
amp = grav*A*np.cos(k*l - frec*(t+dt) + phase)
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + amp

17
src/Mod/Ship/simRun/Sim/initialization.py
View File

@ -40,6 +40,10 @@ class simInitialization:
self.queue = queue
self.loadData(FSmesh, waves)
self.execute()
# Compute time step
self.dt = 0.1
for w in self.waves['data']:
self.dt = np.min(self.dt, w[1]/200.0)
def loadData(self, FSmesh, waves):
""" Convert data to numpy format.
@ -51,8 +55,6 @@ class simInitialization:
nW = len(waves)
# Mesh data
p = np.ndarray((nx,ny, 3), dtype=np.float32)
v = np.ndarray((nx,ny, 3), dtype=np.float32)
f = np.ndarray((nx,ny, 3), dtype=np.float32)
n = np.ndarray((nx,ny, 3), dtype=np.float32)
a = np.ndarray((nx,ny), dtype=np.float32)
phi = np.ndarray((nx,ny), dtype=np.float32)
@ -67,12 +69,6 @@ class simInitialization:
p[i,j,0] = pos.x
p[i,j,1] = pos.y
p[i,j,2] = pos.z
v[i,j,0] = 0.
v[i,j,1] = 0.
v[i,j,2] = 0.
f[i,j,0] = 0.
f[i,j,1] = 0.
f[i,j,2] = 0.
n[i,j,0] = normal.x
n[i,j,1] = normal.y
n[i,j,2] = normal.z
@ -81,9 +77,8 @@ class simInitialization:
Phi[i,j] = 0.
s[i,j] = 0.
ss[i,j] = 0.
self.fs = {'Nx':nx, 'Ny':ny, 'pos':p, 'vel':v, 'acc':f, \
'normal':n, 'area':a, 'velPot':phi, 'accPot':Phi, \
'velSrc':s, 'accSrc':ss}
self.fs = {'Nx':nx, 'Ny':ny, 'pos':p, 'normal':n, 'area':a, \
'velPot':phi, 'accPot':Phi, 'velSrc':s, 'accSrc':ss}
# Waves data
w = np.ndarray((nW, 4), dtype=np.float32)
for i in range(0,nW):

8
src/Mod/Ship/simRun/Simulation.py
View File

@ -87,9 +87,12 @@ class FreeCADShipSimulation(threading.Thread):
init = simInitialization(self.FSmesh,self.waves,self.context,self.queue)
matGen = simMatrixGen(self.context,self.queue)
solver = simComputeSources(self.context,self.queue)
fsEvol = simFSEvolution(self.context,self.queue)
A = init.A
FS = init.fs
waves = init.waves
dt = init.dt
t = 0.0
msg = Translator.translate("\t[Sim]: Iterating...\n")
FreeCAD.Console.PrintMessage(msg)
while self.active:
@ -99,6 +102,11 @@ class FreeCADShipSimulation(threading.Thread):
msg = Translator.translate("\t\t[Sim]: Solving linear systems...\n")
FreeCAD.Console.PrintMessage(msg)
solver.execute(FS, A)
msg = Translator.translate("\t\t[Sim]: Time integrating...\n")
FreeCAD.Console.PrintMessage(msg)
fsEvol.execute(FS, waves, dt, t)
t = t + dt
FreeCAD.Console.PrintMessage('t = %g s' % (t))
# Set thread as stopped (and prepare it to restarting)
self.active = False
threading.Event().set()