diff --git a/src/Mod/Ship/simRun/Sim/fsEvolution.py b/src/Mod/Ship/simRun/Sim/fsEvolution.py index a4d2e0161..14fdf4ba3 100644 --- a/src/Mod/Ship/simRun/Sim/fsEvolution.py +++ b/src/Mod/Ship/simRun/Sim/fsEvolution.py @@ -45,33 +45,19 @@ class simFSEvolution: @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 + # In order to improve results in really long simulations free surface + # will performed considering external waves and second order effects + # in two different ways. First external waves at time t will be + # substracted, then second order waves will be computed, and finally + # external waves at t+dt will be added. for i in range(0,nx): for j in range(0,ny): - # Get value at pos using characteristics method - gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j]) - gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal) - self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*gradVal - # Velocity potential - self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \ - dt*self.fs['accPot'][i,j] + \ - 0.5*dt*dt*grav*self.fs['pos'][i,j][2] - # Acceleration potential. This is really hard to simulate - # accurately due to numerical diffusion of the function, so - # external waves, and diffracted waves will be computed - # in two different ways: - # * External waves will be considered analitically, - # substracting waves at t, and adding waves at t+dt - # * Second order waves will be computed substracting external - # waves to free surface height, and then imposing boundary - # condition. pos = np.copy(self.fs['pos'][i,j]) + # Substract external waves at time t. for w in waves['data']: A = w[0] T = w[1] @@ -81,16 +67,39 @@ class simFSEvolution: k = 2.0*np.pi/wl frec = 2.0*np.pi/T l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading) - # Substract external waves height in order to know second - # order waves free surface amplitude. - amp = A*np.sin(k*l - frec*(t+dt) + phase) + amp = A*np.sin(k*l - frec*t + phase) pos[2] = pos[2] - amp - # Compute analitic external waves acceleration potential - amp0 = grav*A*np.cos(k*l - frec*t + phase) - amp1 = grav*A*np.cos(k*l - frec*(t+dt) + phase) - self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - amp0 + amp1 - # Now impose free surface boundary condition - # self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + grav*pos[2] + amp = - grav/frec*A*np.sin(k*l - frec*t + phase) + self.fs['velPot'][i,j] = self.fs['velPot'][i,j] - amp + amp = grav*A*np.cos(k*l - frec*t + phase) + self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - amp + # Now compute second order waves using position copy, + # where external waves are excluded, in order impose + # free surface boundary condition relative to second + # order phenomena. + self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \ + dt*self.fs['accPot'][i,j] + # self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + \ + # grav*pos[2] + # Restore external waves to velocity and acceleration + # potentials. + 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 + l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading) + 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 + # Update free surface point position + gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j]) + gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal) + self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*gradVal # Impose values at beach (far free surface) for i in range(0,nx): for j in [0,ny-1]: @@ -107,15 +116,12 @@ class simFSEvolution: ny = self.fs['Ny'] nF = nx*ny grad = np.ndarray((nF,3), dtype=np.float32) - FF = open('gradient', 'w') 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) gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j]) gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal) - FF.write('%g\t%g\n' % (pos[1], gradVal)) - FF.close() return grad def gradientphi(self, pos): diff --git a/src/Mod/Ship/simRun/Sim/fsEvolution.pyc b/src/Mod/Ship/simRun/Sim/fsEvolution.pyc deleted file mode 100644 index e8cf2ab57..000000000 Binary files a/src/Mod/Ship/simRun/Sim/fsEvolution.pyc and /dev/null differ diff --git a/src/Mod/Ship/simRun/Simulation.py b/src/Mod/Ship/simRun/Simulation.py index 6ddd1530a..aa8373df8 100644 --- a/src/Mod/Ship/simRun/Simulation.py +++ b/src/Mod/Ship/simRun/Simulation.py @@ -100,7 +100,6 @@ class FreeCADShipSimulation(threading.Thread): ny = FS['Ny'] msg = Translator.translate("\t[Sim]: Iterating...\n") FreeCAD.Console.PrintMessage(msg) - count = 0 while self.active and self.t < self.endTime: msg = Translator.translate("\t\t[Sim]: Generating linear system matrix...\n") FreeCAD.Console.PrintMessage(msg) @@ -113,13 +112,6 @@ class FreeCADShipSimulation(threading.Thread): fsEvol.execute(FS, waves, dt, self.t) self.t = self.t + dt FreeCAD.Console.PrintMessage('t = %g s\n' % (self.t)) - count = count+1 - FF = open('%d' % (count), 'w') - i=1 - for j in range(0,ny): - FF.write("%g\t%g\t%g\t%g\n" % (FS['pos'][i,j,1], FS['pos'][i,j,2], - FS['velPot'][i,j], FS['accPot'][i,j])) - FF.close() # Set thread as stopped (and prepare it to restarting) self.active = False threading.Event().set() diff --git a/src/Mod/Ship/simRun/Simulation.pyc b/src/Mod/Ship/simRun/Simulation.pyc deleted file mode 100644 index bc8596d8d..000000000 Binary files a/src/Mod/Ship/simRun/Simulation.pyc and /dev/null differ