Redaction de la condition de stabilite pour le stretching.

HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py

@@ -44,13 +44,14 @@ class RK2Stretch(ODESolver):
         @param t : current time.
         @param dt : time step.
         """
-        resultTmp = [np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)]
+        resultTmp = np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
+        K1=np.asarray(fd)
+        K1[:,:,:,:]= K1[:,:,:,:] *0.5*dt + field1[:,:,:] 
         if False in [(field1[i][...].shape == field2[i][...].shape) for i in xrange(dimension)] :
             print "Error, not the same mesh on field1 and  on field2"
-        f.discreteOperator = DifferentialOperator_d(opDif=f, field1= ([field1[i][...] + (fd[i][...] * dt * 0.5 )for i in xrange(dimension)]), field2=field2, result=resultTmp, choice='divergenceProduct')
+        f.discreteOperator = DifferentialOperator_d(opDif=f, field1=K1, field2=field2, result=resultTmp, choice='divergenceProduct')
         f.discreteOperator.apply()
-        for i in range(dimension):
-            res_vorticity[i][...] = field1[i][...] + resultTmp[i][...] * dt
+        res_vorticity[:,:,:] = field1[:,:,:] + resultTmp[:,:,:,:] * dt
         return res_vorticity
 
 

HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py

@@ -42,30 +42,25 @@ class RK4Stretch(ODESolver):
         @param t : current time.
         @param dt : time step.
         """
-        resultTmp = [np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)]
-        K1=fd
-        K2=[np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)]
-        K3=[np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)]
-        K4=[np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)]
-        for i in range(dimension):
-            K1[i][...] =K1[i][...]* dt
-            resultTmp[i][...] = K1[i][...] * 0.5 + field1[i][...]
-        f.discreteOperator = DifferentialOperator_d(f, resultTmp, field2, K2, choice='divergenceProduct')
+        resultTmp = np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
+        K1=np.asarray(fd)
+        K2=np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
+        K3=np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
+        K4=np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
+        K1[:,:,:,:] =K1[:,:,:,:]* dt
+        resultTmp[:,:,:,:] = K1[:,:,:,:] * 0.5 + field1[:,:,:]
+        f.discreteOperator = DifferentialOperator_d(f,resultTmp , field2, K2, choice='divergenceProduct')
         f.discreteOperator.apply()
-        for i in range(dimension):
-            K2[i][...] = K2[i][...] * dt
-            resultTmp[i][...] = K2[i][...] * 0.5 + field1[i][...]
+        K2[:,:,:,:] = K2[:,:,:,:] * dt
+        resultTmp[:,:,:,:] = K2[:,:,:,:] * 0.5 + field1[:,:,:]
         f.discreteOperator = DifferentialOperator_d(f, resultTmp, field2, K3, choice='divergenceProduct')
         f.discreteOperator.apply()
-        for i in range(dimension):
-            K3[i][...] = K3[i][...] * dt
-            resultTmp[i][...] = K3[i][...] + field1[i][...]
+        K3[:,:,:,:] = K3[:,:,:,:] * dt
+        resultTmp[:,:,:,:] = K3[:,:,:,:] + field1[:,:,:]
         f.discreteOperator = DifferentialOperator_d(f, resultTmp, field2, K4, choice='divergenceProduct')
         f.discreteOperator.apply()
-        for i in range(dimension):
-            K4[i][...] = K4[i][...] * dt
-        for i in range(dimension):
-            res_vorticity[i][...] = field1[i][...] + 1./6. * (K1[i][...] + 2.0 * K2[i][...] + 2.0  * K3[i][...] + K4[i][...])
+        K4[:,:,:,:] = K4[:,:,:,:] * dt
+        res_vorticity[:,:,:] = field1[:,:,:] + 1./6. * (K1[:,:,:,:] + 2.0 * K2[:,:,:,:] + 2.0  * K3[:,:,:,:] + K4[:,:,:,:])
         return res_vorticity
 
 

HySoP/hysop/test/test_particular_solvers/test_RK.py

@@ -26,21 +26,21 @@ class test_Stretching(unittest.TestCase):
         return vx, vy, vz
 
     def vorticite(self, x, y, z):
-        wx = self.t+1.+ 0*y*z*x
-        wy = self.t+1.+ 0*y*z*x
-        wz = self.t+1.+ 0*y*z*x
+        wx = self.t+1.
+        wy = self.t+1.
+        wz = self.t+1.
         return wx, wy, wz
 
     def analyticalSolution(self, x, y, z):
-        sx = (0.5*(self.t*self.t)+self.t) * np.cos(x) + 0*y*z
-        sy = (0.5*(self.t*self.t)+self.t) * np.cos(y) + 0*x*z
-        sz = (0.5*(self.t*self.t)+self.t) * np.cos(z) + 0*y*x
+        sx = (0.5*(self.t*self.t)+self.t) * np.cos(x)
+        sy = (0.5*(self.t*self.t)+self.t) * np.cos(y)
+        sz = (0.5*(self.t*self.t)+self.t) * np.cos(z)
         return sx, sy, sz
 
     def testOperatorStretching(self):
         # Parameters
-        nb = 128
-        timeStep = 0.005
+        nb = 32
+        timeStep = 0.0025
         finalTime = 0.01
 
         t0 = time.time()
@@ -66,23 +66,21 @@ class test_Stretching(unittest.TestCase):
 
         ## Setting solver to Problem
         pb.setSolver(finalTime, timeStep, solver_type='basic')
-        Analvorty=[np.zeros([nb, nb, nb], dtype=PARMES_REAL, order=ORDER) for d in xrange(3)]
-        for i in xrange(3):
-            Analvorty[i]= self.vorticite(topo3D.mesh.coords[0], \
+        Analvorty=np.asarray([np.zeros([nb, nb, nb], dtype=PARMES_REAL, order=ORDER) for d in xrange(3)])
+        Analvorty[:,:,:,:]= np.vectorize(self.vorticite)(topo3D.mesh.coords[0], \
                                           topo3D.mesh.coords[1], \
-                                          topo3D.mesh.coords[2])[i]
+                                          topo3D.mesh.coords[2])[:]
 
-        pb.solver.ODESolver= RK4Stretch #  RK2Stretch # Euler#
+        pb.solver.ODESolver=   RK2Stretch # Euler#RK4Stretch #
         pb.initSolver()
         t1 = time.time()
 
         ## Solve problem
         pb.solve()
         self.t = pb.timer.t
-        for i in range(3):
-            Analvorty[i]=Analvorty[i]+self.analyticalSolution(topo3D.mesh.coords[0], \
+        Analvorty[:,:,:,:]=Analvorty[:,:,:,:]+np.vectorize(self.analyticalSolution)(topo3D.mesh.coords[0], \
                                             topo3D.mesh.coords[1], \
-                                            topo3D.mesh.coords[2])[i]
+                                            topo3D.mesh.coords[2])[:]
 
 ## Visual comparison between the numerical resolution and the analytical resolution
         for i in range(1):
@@ -95,8 +93,8 @@ class test_Stretching(unittest.TestCase):
             plt.axis([-2*np.pi,2*np.pi,0,max(abs(Analvorty[i][:,0,0] - vorti.discreteField[i].data[i][:,0,0]))])
             plt.plot(topo3D.mesh.coords[0][:,0,0], abs(Analvorty[i][:,0,0] - vorti.discreteField[i].data[i][:,0,0]), '-' )
             plt.legend([u"Erreur"])
-#            plt.show()
-            npt.assert_array_less(abs(Analvorty[i][...] - vorti.discreteField[i].data[i][...]),0.002)
+            plt.show()
+        npt.assert_array_less(abs(Analvorty[:,:,:,:] - vorti.discreteField[0].data[:,:,:]),0.002)
         tf = time.time()
 
         print "\n"