Improved free surface evolving method in order to support long simulations.

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):