diff --git a/src/Mod/Ship/CMakeLists.txt b/src/Mod/Ship/CMakeLists.txt index 68fc8e924..300f999f8 100644 --- a/src/Mod/Ship/CMakeLists.txt +++ b/src/Mod/Ship/CMakeLists.txt @@ -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}) diff --git a/src/Mod/Ship/Makefile.am b/src/Mod/Ship/Makefile.am index 3a8c9718d..308131bba 100644 --- a/src/Mod/Ship/Makefile.am +++ b/src/Mod/Ship/Makefile.am @@ -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) diff --git a/src/Mod/Ship/simRun/Sim/__init__.py b/src/Mod/Ship/simRun/Sim/__init__.py index a2fd6db39..aabf4a621 100644 --- a/src/Mod/Ship/simRun/Sim/__init__.py +++ b/src/Mod/Ship/simRun/Sim/__init__.py @@ -24,3 +24,4 @@ from initialization import * from matrixGen import * from computeSources import * +from fsEvolution import * diff --git a/src/Mod/Ship/simRun/Sim/fsEvolution.py b/src/Mod/Ship/simRun/Sim/fsEvolution.py new file mode 100644 index 000000000..cd68fe085 --- /dev/null +++ b/src/Mod/Ship/simRun/Sim/fsEvolution.py @@ -0,0 +1,136 @@ +#*************************************************************************** +#* * +#* Copyright (c) 2011, 2012 * +#* Jose Luis Cercos Pita * +#* * +#* 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 + \ No newline at end of file diff --git a/src/Mod/Ship/simRun/Sim/initialization.py b/src/Mod/Ship/simRun/Sim/initialization.py index 823ff535e..b2da085ff 100644 --- a/src/Mod/Ship/simRun/Sim/initialization.py +++ b/src/Mod/Ship/simRun/Sim/initialization.py @@ -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): diff --git a/src/Mod/Ship/simRun/Simulation.py b/src/Mod/Ship/simRun/Simulation.py index a0bb2ff37..598b6178e 100644 --- a/src/Mod/Ship/simRun/Simulation.py +++ b/src/Mod/Ship/simRun/Simulation.py @@ -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()