diff --git a/Docs/Stretching.pdf b/Docs/Stretching.pdf
index 25b8fe0be4fba2f4276e345e7145982c370428ba..34d131b52353ef324b687bb0651efad887b63a6e 100644
Binary files a/Docs/Stretching.pdf and b/Docs/Stretching.pdf differ
diff --git a/Docs/Stretching.tex b/Docs/Stretching.tex
index 6e3226adc57542cb1cba4cd763f05fd2e99b5f10..ac88217e3ad55fa5647b11bb44fc64d1ca38eed1 100644
--- a/Docs/Stretching.tex
+++ b/Docs/Stretching.tex
@@ -240,114 +240,7 @@ Adams-Moulton 3 & 0.44984132 & 3 \\
-\section{Validation Sch\'ema Runge-Kutta}
-Pour valider les sch\'emas de Runge-Kutta, nous sommes partis de l'exemple suivant~:
-\begin{align}
- \overrightarrow{u} =
-\left(\begin{array}{l}
-\sin(x)\\
-\sin(y)\\
-\sin(z)
-\end{array}\right)&&\\
- \overrightarrow{w} =
-\left(\begin{array}{l}
-1+t\\
-1+t\\
-1+t
-\end{array}\right)&&
-\end{align}
-
-La solution analytique est alors~:
-\begin{equation}
- \omega_\text{analytique} =
-\left(\begin{array}{l}
-(\frac{t^2}{2} + t) \cos(x)\\
-(\frac{t^2}{2} + t) \cos(y)\\
-(\frac{t^2}{2} + t) \cos(z)
-\end{array}\right)
-\end{equation}
-
-Seulement cette solution analytique n'est valable que pour le sch\'ema d'Euler. Pour les sch\'emas de type Runge-Kutta, il faut rajouter un terme correcteur d\^u aux erreurs de d\'erivation aux diff\'erents temps.\\
-On a alors~:
-\begin{equation}
- \omega = \omega_\text{analytique} + F(t)
-\end{equation}
-
-Les fonctions correctrices pour $t^1$ pour cet exemple sont les suivantes~:
-\begin{tabular}{|c|c|}
-\hline
-Sch\'ema & Terme correcteur\\
- \hline
-RK 2 & $\left(\begin{array}{l}
-\cos^2(x)(\frac{\Delta t^2}{2} + \frac{\Delta t^3}{4})\\
-\cos^2(y)(\frac{\Delta t^2}{2} + \frac{\Delta t^3}{4})\\
-\cos^2(z)(\frac{\Delta t^2}{2} + \frac{\Delta t^3}{4})
-\end{array}\right)
-$\\ \hline
-RK 3 & $\left(\begin{array}{l}
-\cos^2(x)(\frac{\Delta t^2}{2} + \frac{\Delta t^3}{2})+\cos^3(x)(\frac{\Delta t^3}{6} + \frac{\Delta t^4}{6})\\
-\cos^2(y)(\frac{\Delta t^2}{2} + \frac{\Delta t^3}{2})+\cos^3(y)(\frac{\Delta t^3}{6} + \frac{\Delta t^4}{6})\\
-\cos^2(z)(\frac{\Delta t^2}{2} + \frac{\Delta t^3}{2})+\cos^3(z)(\frac{\Delta t^3}{6} + \frac{\Delta t^4}{6})
-\end{array}\right)
-$\\ \hline
-RK 4 & $\left(\begin{array}{l}
-\frac{1}{6} \cos(x)\biggl(-\Delta t^2 + \cos(x)\Delta t^2(3 + \frac{7}{2}\Delta t)\\ + \frac{\Delta t^3}{2}\cos^2(x)(2+3\Delta t) + \frac{\Delta t^4}{4} \cos^3(x)(1 + 2\Delta t)\biggr)\\
-\frac{1}{6} \cos(y)\biggl(-\Delta t^2 + \cos(y)\Delta t^2(3 + \frac{7}{2}\Delta t)\\ + \frac{\Delta t^3}{2}\cos^2(y)(2+3\Delta t) + \frac{\Delta t^4}{4} \cos^3(y)(1 + 2\Delta t)\biggr)\\
-\frac{1}{6} \cos(z)\biggl(-\Delta t^2 + \cos(z)\Delta t^2(3 + \frac{7}{2}\Delta t) \\+ \frac{\Delta t^3}{2}\cos^2(z)(2+3\Delta t) + \frac{\Delta t^4}{4} \cos^3(z)(1 + 2\Delta t)\biggr)
-\end{array}\right)
-$\\ \hline
-\end{tabular}
-Avec les termes correcteurs nous obtenons les erreurs suivantes figures(\ref{fig:errorRK2}-\ref{fig:errorRK3}-\ref{fig:errorRK3TVD}-\ref{fig:errorRK4}).
-\begin{figure}[ht]
- \centering
-\includegraphics[width=1\linewidth]{RK2_error}
-\caption{\label{fig:errorRK2} \`A gauche la diff\'erencre entre la solution num\'erique pour le sch\'ema de Runge-Kutta 2 et analytique sans terme correctif, \`a droite avec terme correctif.}
-\end{figure}
-
-\begin{figure}[ht]
- \centering
-\includegraphics[width=1\linewidth]{RK3_error}
-\caption{\label{fig:errorRK3} \`A gauche la diff\'erencre entre la solution num\'erique pour le sch\'ema de Runge-Kutta 3 et analytique sans terme correctif, \`a droite avec terme correctif.}
-\end{figure}
-\begin{figure}[ht]
- \centering
-\includegraphics[width=1\linewidth]{RK3_TVD_error}
-\caption{\label{fig:errorRK3TVD} \`A gauche la diff\'erencre entre la solution num\'erique pour le sch\'ema de Runge-Kutta 3 TVD et analytique sans terme correctif, \`a droite avec terme correctif.}
-\end{figure}
-\begin{figure}[ht]
- \centering
-\includegraphics[width=1\linewidth]{RK4_error}
-\caption{\label{fig:errorRK4} \`A gauche la diff\'erencre entre la solution num\'erique pour le sch\'ema de Runge-Kutta 4 et analytique sans terme correctif, \`a droite avec terme correctif.}
-\end{figure}
-
-Les diff\'erents sch\'emas de Runge-Kutta 3 que l'on a impl\'ement\'e~:
-\begin{itemize}
- \item[Version 1]
-$\begin{array}{ll}
-\omega^{n+1} = \omega^n + \frac{1}{4}(K_1 + 3K_3)\\
-K_1= \Delta t f(u^n,\omega^n,t^n)\\
-K_2= \Delta t f(u^n,\omega^n + \frac{K_1}{3},t^n + \frac{1}{3}\Delta t)\\
-K_3= \Delta t f(u^n,\omega^n+ \frac{2K_2}{3},t^n+ \frac{2}{3}\Delta t)
-\end{array}
-$
- \item[Version 2]
-$\begin{array}{ll}
-\omega^{n+1} = \omega^n + \frac{1}{6}(K_1 + 4K_2+ K_3)\\
-K_1= \Delta t f(u^n,\omega^n,t^n)\\
-K_2= \Delta t f(u^n,\omega^n + \frac{K_1}{2},t^n + \frac{1}{2}\Delta t)\\
-K_3= \Delta t f(u^n,\omega^n +2K_2 -K_1,t^n + \Delta t)
-\end{array}
-$
- \item[Version TVD]
-$\begin{array}{ll}
-\omega^{n+1} = \omega^n + \frac{2}{3}(K_2 + K_3)\\
-K_1= \Delta t f(u^n,\omega^n,t^n)\\
-K_2= \Delta t f(u^n,\frac{3}{4}\omega^n + \frac{K_1}{4},t^n + \Delta t)\\
-K_3= \Delta t f(u^n,\omega^n + K_2,t^n + \frac{1}{2}\Delta t)
-\end{array}
-$
-\end{itemize}
\end{document}
diff --git a/Docs/images/RK2_error.png b/Docs/images/RK2_error.png
deleted file mode 100644
index 6078314152230b96b90fb0495c7958953d3ebcf7..0000000000000000000000000000000000000000
Binary files a/Docs/images/RK2_error.png and /dev/null differ
diff --git a/Docs/images/RK3_TVD_error.png b/Docs/images/RK3_TVD_error.png
deleted file mode 100644
index 7c1cb870b3d8a4b8434bad947a5aa0cef5ec4edf..0000000000000000000000000000000000000000
Binary files a/Docs/images/RK3_TVD_error.png and /dev/null differ
diff --git a/Docs/images/RK3_error.png b/Docs/images/RK3_error.png
deleted file mode 100644
index 21400693693424e40a3d307832bd923dc7ef839b..0000000000000000000000000000000000000000
Binary files a/Docs/images/RK3_error.png and /dev/null differ
diff --git a/Docs/images/RK4_error.png b/Docs/images/RK4_error.png
deleted file mode 100644
index 4076b62bce759b10d45446ac82cfb74f705619cc..0000000000000000000000000000000000000000
Binary files a/Docs/images/RK4_error.png and /dev/null differ
diff --git a/HySoP/hysop/operator/differentialOperator.py b/HySoP/hysop/operator/differentialOperator.py
index e4337174417c269af3265a2f088c67560e6a9a87..46888ddefef1b33645e70869ab25fc4885ebe0de 100755
--- a/HySoP/hysop/operator/differentialOperator.py
+++ b/HySoP/hysop/operator/differentialOperator.py
@@ -13,13 +13,14 @@ class DifferentialOperator(ContinuousOperator):
Differential operator
"""
- def __init__(self, field1, field2, choice=None): # changer commentaires
+ def __init__(self, field1, field2, choice=None):
"""
Constructor.
Create a Differential operator from given velocity and vorticity variables.
- @param velocity ContinuousVectorField : velocity variable.
- @param vorticity ContinuousVectorField : vorticity variable.
+ @param field1 : velocity variable.
+ @param field2 : vorticity variable.
+ @param choice : choice of the differential operation in : curl , divWU , gradU
"""
ContinuousOperator.__init__(self)
self.field1 = field1
@@ -27,19 +28,18 @@ class DifferentialOperator(ContinuousOperator):
self.choice = choice
self.addVariable([field1, field2])
- def discretize(self, field1=None, field2=None, result=None, maxgersh=None): # changer commentaires
+ def discretize(self, field1=None, field2=None, resolution=None, meshSize=None):
"""
Stretching operator discretization method.
Create a discrete Stretching operator from given specifications.
- @param idVelocityD : Index of velocity discretisation to use.
- @param idVorticityD : Index of vorticity discretisation to use.
- @param res_vorticity : result stretch (vorticity).
- @param method : the method to use.
- @param maxgersh : result Gershgorin's condition.
+ @param field1 : Index of velocity discretisation to use.
+ @param field2 : Index of vorticity discretisation to use.
+ @param resolution : local resolution.
+ @param meshSize : mesh size.
"""
if self.discreteOperator is None:
- self.discreteOperator = DifferentialOperator_d(self, field1, field2, result, self.choice, maxgersh)
+ self.discreteOperator = DifferentialOperator_d(field1, field2, self.choice, resolution, meshSize)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/differentialOperator_d.py b/HySoP/hysop/operator/differentialOperator_d.py
index 3200094aad4c8d94de71f6ab2a3a6bbab82d7da7..247ad547285eeaef326cfcfd53cd3f980c919586 100755
--- a/HySoP/hysop/operator/differentialOperator_d.py
+++ b/HySoP/hysop/operator/differentialOperator_d.py
@@ -15,7 +15,7 @@ class DifferentialOperator_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, opDif, field1, field2, result=None, choice=None, maxgersh=None):
+ def __init__(self, field1, field2, choice=None, resolution=None, meshSize=None):
"""
Constructor.
Create a Stretching operator on a given continuous domain.
@@ -26,295 +26,180 @@ class DifferentialOperator_d(DiscreteOperator):
DiscreteOperator.__init__(self)
self.field1 = field1
self.field2 = field2
- self.result = result
- self.maxgersh = maxgersh
self.choice = choice
+ self.resolution = resolution
+ self.meshSize = meshSize
self.compute_time = 0.
def apply(self):
"""
Apply operator.
"""
- if False in [(self.field1[i][...].shape == self.field2[i][...].shape) for i in xrange(self.field2.topology.dim)] :
- print "Error, not the same topology for field1 and field2"
- sys.exit(0)
- else:
- mesh_size = self.field2.topology.mesh.size
- resolution = self.field2.topology.mesh.resolution
- ind=np.arange(resolution[0])
-
- if self.choice == 'masseConservation':
+# if False in [(self.field1[i][...].shape == self.field2[i][...].shape) for i in xrange(self.field2.topology.dim)] :
+# print "Error, not the same topology for field1 and field2"
+# sys.exit(0)
+# else:
+# mesh_size = self.field2.topology.mesh.size
+# resolution = self.field2.topology.mesh.resolution
+ ind=np.arange(self.resolution[0])
+
+ if self.choice == 'divWU':
## second order scheme
## X components of temp and result
# temp1 = ( self.field1[0][np.roll(ind,-1),...] * self.field2[0][np.roll(ind,-1),...] -\
-# self.field1[0][np.roll(ind,+1),...] * self.field2[0][np.roll(ind,+1),...]) / (2. * mesh_size[0])
+# self.field1[0][np.roll(ind,+1),...] * self.field2[0][np.roll(ind,+1),...]) / (2. * self.meshSize[0])
# temp2 = (self.field1[1][:,np.roll(ind,-1),:] * self.field2[0][:,np.roll(ind,-1),:] -\
-# self.field1[1][:,np.roll(ind,+1),:] * self.field2[0][:,np.roll(ind,+1),:] ) / (2. * mesh_size[1])
+# self.field1[1][:,np.roll(ind,+1),:] * self.field2[0][:,np.roll(ind,+1),:] ) / (2. * self.meshSize[1])
#
# temp3 = (self.field1[2][...,np.roll(ind,-1)] * self.field2[0][...,np.roll(ind,-1)] -\
-# self.field1[2][...,np.roll(ind,+1)] * self.field2[0][...,np.roll(ind,+1)] ) / (2. * mesh_size[2])
+# self.field1[2][...,np.roll(ind,+1)] * self.field2[0][...,np.roll(ind,+1)] ) / (2. * self.meshSize[2])
# self.result[0][...]= temp1 + temp2 + temp3
## Y components of temp and result
# temp1 = (self.field1[0][np.roll(ind,-1),...] * self.field2[1][np.roll(ind,-1),...] -\
-# self.field1[0][np.roll(ind,+1),...] * self.field2[1][np.roll(ind,+1),...] ) / (2. * mesh_size[0])
+# self.field1[0][np.roll(ind,+1),...] * self.field2[1][np.roll(ind,+1),...] ) / (2. * self.meshSize[0])
# temp2 = (self.field1[1][:,np.roll(ind,-1),:] * self.field2[1][:,np.roll(ind,-1),:] -\
-# self.field1[1][:,np.roll(ind,+1),:] * self.field2[1][:,np.roll(ind,+1),:] ) / (2. * mesh_size[1])
+# self.field1[1][:,np.roll(ind,+1),:] * self.field2[1][:,np.roll(ind,+1),:] ) / (2. * self.meshSize[1])
#
# temp3 = (self.field1[2][...,np.roll(ind,-1)] * self.field2[1][...,np.roll(ind,-1)] -\
-# self.field1[2][...,np.roll(ind,+1)] * self.field2[1][...,np.roll(ind,+1)] ) / (2. * mesh_size[2])
+# self.field1[2][...,np.roll(ind,+1)] * self.field2[1][...,np.roll(ind,+1)] ) / (2. * self.meshSize[2])
# self.result[1][...]= temp1 + temp2 +temp3
## Z components of temp and result
# temp1 = (self.field1[0][np.roll(ind,-1),...] * self.field2[2][np.roll(ind,-1),...] -\
-# self.field1[0][np.roll(ind,+1),...] * self.field2[2][np.roll(ind,+1),...] ) / (2. * mesh_size[0])
+# self.field1[0][np.roll(ind,+1),...] * self.field2[2][np.roll(ind,+1),...] ) / (2. * self.meshSize[0])
# temp2 = (self.field1[1][:,np.roll(ind,-1),:] * self.field2[2][:,np.roll(ind,-1),:] -\
-# self.field1[1][:,np.roll(ind,+1),:] * self.field2[2][:,np.roll(ind,+1),:] ) / (2. * mesh_size[1])
+# self.field1[1][:,np.roll(ind,+1),:] * self.field2[2][:,np.roll(ind,+1),:] ) / (2. * self.meshSize[1])
#
# temp3 = (self.field1[2][...,np.roll(ind,-1)] * self.field2[2][...,np.roll(ind,-1)] -\
-# self.field1[2][...,np.roll(ind,+1)] * self.field2[2][...,np.roll(ind,+1)]) / (2. * mesh_size[2])
+# self.field1[2][...,np.roll(ind,+1)] * self.field2[2][...,np.roll(ind,+1)]) / (2. * self.meshSize[2])
# X components of temp and result
temp1 = (1.0 * self.field1[0][np.roll(ind,+2),...] * self.field2[0][np.roll(ind,+2),...] -
8.0 * self.field1[0][np.roll(ind,+1),...] * self.field2[0][np.roll(ind,+1),...] +\
8.0 * self.field1[0][np.roll(ind,-1),...] * self.field2[0][np.roll(ind,-1),...] -\
- 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[0][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+ 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[0][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
temp2 = (1.0 * self.field1[1][:,np.roll(ind,+2),:] * self.field2[0][:,np.roll(ind,+2),:] -\
8.0 * self.field1[1][:,np.roll(ind,+1),:] * self.field2[0][:,np.roll(ind,+1),:] +\
8.0 * self.field1[1][:,np.roll(ind,-1),:] * self.field2[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[0][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+ 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[0][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
temp3 = (1.0 * self.field1[2][...,np.roll(ind,+2)] * self.field2[0][...,np.roll(ind,+2)] -\
8.0 * self.field1[2][...,np.roll(ind,+1)] * self.field2[0][...,np.roll(ind,+1)] +\
8.0 * self.field1[2][...,np.roll(ind,-1)] * self.field2[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[0][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+ 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[0][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
- self.result[0][...]= temp1 + temp2 + temp3
+ tmp1 = np.array([temp1 + temp2 + temp3])
# Y components of temp and result
temp1 = (1.0 * self.field1[0][np.roll(ind,+2),...] * self.field2[1][np.roll(ind,+2),...] -
8.0 * self.field1[0][np.roll(ind,+1),...] * self.field2[1][np.roll(ind,+1),...] +\
8.0 * self.field1[0][np.roll(ind,-1),...] * self.field2[1][np.roll(ind,-1),...] -\
- 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[1][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+ 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[1][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
temp2 = (1.0 * self.field1[1][:,np.roll(ind,+2),:] * self.field2[1][:,np.roll(ind,+2),:] -\
8.0 * self.field1[1][:,np.roll(ind,+1),:] * self.field2[1][:,np.roll(ind,+1),:] +\
8.0 * self.field1[1][:,np.roll(ind,-1),:] * self.field2[1][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[1][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+ 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[1][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
temp3 = (1.0 * self.field1[2][...,np.roll(ind,+2)] * self.field2[1][...,np.roll(ind,+2)] -\
8.0 * self.field1[2][...,np.roll(ind,+1)] * self.field2[1][...,np.roll(ind,+1)] +\
8.0 * self.field1[2][...,np.roll(ind,-1)] * self.field2[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[1][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+ 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[1][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
- self.result[1][...]= temp1 + temp2 +temp3
+ tmp2 = np.array([temp1 + temp2 + temp3])
# Z components of temp and result
temp1 = (1.0 * self.field1[0][np.roll(ind,+2),...] * self.field2[2][np.roll(ind,+2),...] -
8.0 * self.field1[0][np.roll(ind,+1),...] * self.field2[2][np.roll(ind,+1),...] +\
8.0 * self.field1[0][np.roll(ind,-1),...] * self.field2[2][np.roll(ind,-1),...] -\
- 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[2][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+ 1.0 * self.field1[0][np.roll(ind,-2),...] * self.field2[2][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
temp2 = (1.0 * self.field1[1][:,np.roll(ind,+2),:] * self.field2[2][:,np.roll(ind,+2),:] -\
8.0 * self.field1[1][:,np.roll(ind,+1),:] * self.field2[2][:,np.roll(ind,+1),:] +\
8.0 * self.field1[1][:,np.roll(ind,-1),:] * self.field2[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[2][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+ 1.0 * self.field1[1][:,np.roll(ind,-2),:] * self.field2[2][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
temp3 = (1.0 * self.field1[2][...,np.roll(ind,+2)] * self.field2[2][...,np.roll(ind,+2)] -\
8.0 * self.field1[2][...,np.roll(ind,+1)] * self.field2[2][...,np.roll(ind,+1)] +\
8.0 * self.field1[2][...,np.roll(ind,-1)] * self.field2[2][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[2][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
-
- self.result[2][...]= temp1 + temp2 +temp3
-
-
- elif self.choice == 'masseConservationF':
-
-# X components of temp and result
- temp1 = (-25./12. * self.field1[0][np.roll(ind,0),...] * self.field2[0][np.roll(ind,0),...] +
- 4. * self.field1[0][np.roll(ind,+1),...] * self.field2[0][np.roll(ind,+1),...] +\
- -3. * self.field1[0][np.roll(ind,+2),...] * self.field2[0][np.roll(ind,+2),...] +\
- 4./3. * self.field1[0][np.roll(ind,+3),...] * self.field2[0][np.roll(ind,+3),...] +\
- -1./4. * self.field1[0][np.roll(ind,+4),...] * self.field2[0][np.roll(ind,+4),...]) / (mesh_size[0])
-
- temp2 = ( -25./12. * self.field1[1][:,np.roll(ind,+0),:] * self.field2[0][:,np.roll(ind,0),:] +\
- 4. * self.field1[1][:,np.roll(ind,+1),:] * self.field2[0][:,np.roll(ind,+1),:] +\
- -3. * self.field1[1][:,np.roll(ind,+2),:] * self.field2[0][:,np.roll(ind,+2),:] +\
- 4./3. * self.field1[1][:,np.roll(ind,+3),:] * self.field2[0][:,np.roll(ind,+3),:] +\
- -1./4. * self.field1[1][:,np.roll(ind,+4),:] * self.field2[0][:,np.roll(ind,+4),:]) / ( mesh_size[1])
-
- temp3 = ( -25./12. * self.field1[2][...,np.roll(ind,+0)] * self.field2[0][...,np.roll(ind,+0)] +\
- 4. * self.field1[2][...,np.roll(ind,+1)] * self.field2[0][...,np.roll(ind,+1)] +\
- -3. * self.field1[2][...,np.roll(ind,+2)] * self.field2[1][...,np.roll(ind,+2)] +\
- 4./3. * self.field1[2][...,np.roll(ind,+3)] * self.field2[0][...,np.roll(ind,+3)] +\
- -1./4. * self.field1[2][...,np.roll(ind,+4)] * self.field2[0][...,np.roll(ind,+4)]) / ( mesh_size[2])
-
- self.result[0][...]= temp1 + temp2 + temp3
-
-# Y components of temp and result
- temp1 = (-25./12. * self.field1[0][np.roll(ind,+0),...] * self.field2[1][np.roll(ind,+0),...] +
- 4. * self.field1[0][np.roll(ind,+1),...] * self.field2[1][np.roll(ind,+1),...] +\
- -3. * self.field1[0][np.roll(ind,+2),...] * self.field2[1][np.roll(ind,+2),...] +\
- 4./3. * self.field1[0][np.roll(ind,+3),...] * self.field2[1][np.roll(ind,+3),...] +\
- -1./4. * self.field1[0][np.roll(ind,+4),...] * self.field2[1][np.roll(ind,+4),...]) / ( mesh_size[0])
-
- temp2 = (-25./12. * self.field1[1][:,np.roll(ind,+0),:] * self.field2[1][:,np.roll(ind,+0),:] +\
- 4. * self.field1[1][:,np.roll(ind,+1),:] * self.field2[1][:,np.roll(ind,+1),:] +\
- -3. * self.field1[1][:,np.roll(ind,+2),:] * self.field2[1][:,np.roll(ind,+2),:] +\
- 4./3. * self.field1[1][:,np.roll(ind,+3),:] * self.field2[1][:,np.roll(ind,+3),:] +\
- -1./4. * self.field1[1][:,np.roll(ind,+4),:] * self.field2[1][:,np.roll(ind,+4),:]) / ( mesh_size[1])
-
- temp3 = (-25./12. * self.field1[2][...,np.roll(ind,+0)] * self.field2[1][...,np.roll(ind,+0)] +\
- 4. * self.field1[2][...,np.roll(ind,+1)] * self.field2[1][...,np.roll(ind,+1)] +\
- -3. * self.field1[2][...,np.roll(ind,+2)] * self.field2[1][...,np.roll(ind,+2)] +\
- 4./3. * self.field1[2][...,np.roll(ind,+3)] * self.field2[1][...,np.roll(ind,+3)] +\
- -1./4. * self.field1[2][...,np.roll(ind,+4)] * self.field2[1][...,np.roll(ind,+4)]) / ( mesh_size[2])
-
- self.result[1][...]= temp1 + temp2 +temp3
-
-# Z components of temp and result
- temp1 = (-25./12. * self.field1[0][np.roll(ind,+0),...] * self.field2[2][np.roll(ind,+0),...] +
- 4. * self.field1[0][np.roll(ind,+1),...] * self.field2[2][np.roll(ind,+1),...] +\
- -3. * self.field1[0][np.roll(ind,+2),...] * self.field2[2][np.roll(ind,+2),...] +\
- 4./3. * self.field1[0][np.roll(ind,+3),...] * self.field2[2][np.roll(ind,+3),...] +\
- -1./4. * self.field1[0][np.roll(ind,+4),...] * self.field2[2][np.roll(ind,+4),...]) / ( mesh_size[0])
-
- temp2 = (-25./12. * self.field1[1][:,np.roll(ind,+0),:] * self.field2[2][:,np.roll(ind,+0),:] +\
- 4. * self.field1[1][:,np.roll(ind,+1),:] * self.field2[2][:,np.roll(ind,+1),:] +\
- -3. * self.field1[1][:,np.roll(ind,+2),:] * self.field2[2][:,np.roll(ind,+2),:] +\
- 4./3. * self.field1[1][:,np.roll(ind,+3),:] * self.field2[2][:,np.roll(ind,+3),:] +\
- -1./4. * self.field1[1][:,np.roll(ind,+4),:] * self.field2[2][:,np.roll(ind,+4),:]) / ( mesh_size[1])
-
- temp3 = (-25./12. * self.field1[2][...,np.roll(ind,+0)] * self.field2[2][...,np.roll(ind,+0)] +\
- 4. * self.field1[2][...,np.roll(ind,+1)] * self.field2[2][...,np.roll(ind,+1)] +\
- -3. * self.field1[2][...,np.roll(ind,+2)] * self.field2[2][...,np.roll(ind,+2)] +\
- 4./3. * self.field1[2][...,np.roll(ind,+3)] * self.field2[2][...,np.roll(ind,+3)] +\
- -1./4. * self.field1[2][...,np.roll(ind,+4)] * self.field2[2][...,np.roll(ind,+4)]) / ( mesh_size[2])
+ 1.0 * self.field1[2][...,np.roll(ind,-2)] * self.field2[2][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
- self.result[2][...]= temp1 + temp2 +temp3
+ tmp3 = np.array([temp1 + temp2 + temp3])
+ result = np.concatenate((np.array([tmp1]), np.array([tmp2]), np.array([tmp3])))
+ return result
- elif self.choice == 'divConservation':
+ elif self.choice == 'gradU':
-## fourth order scheme
+ maxgersh = np.zeros(2, dtype=PARMES_REAL, order=ORDER)
+## Fourth order scheme
# X components of temp and result
temp1 = (1.0 * self.field2[0][np.roll(ind,+2),...] -
8.0 * self.field2[0][np.roll(ind,+1),...] +\
8.0 * self.field2[0][np.roll(ind,-1),...] -\
- 1.0 * self.field2[0][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+ 1.0 * self.field2[0][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
temp2 = (1.0 * self.field2[0][:,np.roll(ind,+2),:] -\
8.0 * self.field2[0][:,np.roll(ind,+1),:] +\
8.0 * self.field2[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[0][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+ 1.0 * self.field2[0][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
temp3 = (1.0 * self.field2[0][...,np.roll(ind,+2)] -\
8.0 * self.field2[0][...,np.roll(ind,+1)] +\
8.0 * self.field2[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[0][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+ 1.0 * self.field2[0][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
- max1= np.max(abs(temp1)+abs(temp2)+abs(temp3))
+ maxstr1= np.max(abs(temp1)+abs(temp2)+abs(temp3))
maxadv1= np.max(abs(temp1))
- self.result[0][...]= temp1 * self.field1[0][...] + temp2 * self.field1[1][...] + temp3 * self.field1[2][...]
# Y components of temp and result
- temp1 = (1.0 * self.field2[1][np.roll(ind,+2),...] -
+ temp4 = (1.0 * self.field2[1][np.roll(ind,+2),...] -
8.0 * self.field2[1][np.roll(ind,+1),...] +\
8.0 * self.field2[1][np.roll(ind,-1),...] -\
- 1.0 * self.field2[1][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+ 1.0 * self.field2[1][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field2[1][:,np.roll(ind,+2),:] -\
+ temp5 = (1.0 * self.field2[1][:,np.roll(ind,+2),:] -\
8.0 * self.field2[1][:,np.roll(ind,+1),:] +\
8.0 * self.field2[1][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[1][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+ 1.0 * self.field2[1][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field2[1][...,np.roll(ind,+2)] -\
+ temp6 = (1.0 * self.field2[1][...,np.roll(ind,+2)] -\
8.0 * self.field2[1][...,np.roll(ind,+1)] +\
8.0 * self.field2[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[1][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+ 1.0 * self.field2[1][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
- max2= np.max(abs(temp1)+abs(temp2)+abs(temp3))
+ maxstr2= np.max(abs(temp1)+abs(temp2)+abs(temp3))
maxadv2= np.max(abs(temp2))
- self.result[1][...]= temp1 * self.field1[0][...] + temp2 * self.field1[1][...] + temp3 * self.field1[2][...]
# Z components of temp and result
- temp1 = (1.0 * self.field2[2][np.roll(ind,+2),...] -
+ temp7 = (1.0 * self.field2[2][np.roll(ind,+2),...] -
8.0 * self.field2[2][np.roll(ind,+1),...] +\
8.0 * self.field2[2][np.roll(ind,-1),...] -\
- 1.0 * self.field2[2][np.roll(ind,-2),...]) / (12. * mesh_size[0])
+ 1.0 * self.field2[2][np.roll(ind,-2),...]) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field2[2][:,np.roll(ind,+2),:] -\
+ temp8 = (1.0 * self.field2[2][:,np.roll(ind,+2),:] -\
8.0 * self.field2[2][:,np.roll(ind,+1),:] +\
8.0 * self.field2[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[2][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
+ 1.0 * self.field2[2][:,np.roll(ind,-2),:]) / (12. * self.meshSize[1])
- temp3 = (1.0 * self.field2[2][...,np.roll(ind,+2)] -\
+ temp9 = (1.0 * self.field2[2][...,np.roll(ind,+2)] -\
8.0 * self.field2[2][...,np.roll(ind,+1)] +\
8.0 * self.field2[2][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[2][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+ 1.0 * self.field2[2][...,np.roll(ind,-2)]) / (12. * self.meshSize[2])
- max3= np.max(abs(temp1)+abs(temp2)+abs(temp3))
+ maxstr3= np.max(abs(temp1)+abs(temp2)+abs(temp3))
maxadv3= np.max(abs(temp3))
- self.result[2][...]= temp1 * self.field1[0][...] + temp2 * self.field1[1][...] + temp3 * self.field1[2][...]
- self.maxgersh[0] = max(max1,max2,max3)
-# self.maxgersh[0] = max(maxadv1,maxadv2,maxadv3)
-
- elif self.choice == 'divConservationGradU':
-
-## fourth order scheme
-# First line of Grad U matrix
- self.result[0][...] = (1.0 * self.field2[0][np.roll(ind,+2),...] -
- 8.0 * self.field2[0][np.roll(ind,+1),...] +\
- 8.0 * self.field2[0][np.roll(ind,-1),...] -\
- 1.0 * self.field2[0][np.roll(ind,-2),...]) / (12. * mesh_size[0])
-
- self.result[1][...] = (1.0 * self.field2[0][:,np.roll(ind,+2),:] -\
- 8.0 * self.field2[0][:,np.roll(ind,+1),:] +\
- 8.0 * self.field2[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[0][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
-
- self.result[2][...] = (1.0 * self.field2[0][...,np.roll(ind,+2)] -\
- 8.0 * self.field2[0][...,np.roll(ind,+1)] +\
- 8.0 * self.field2[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[0][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
-
-# Second line of Grad U matrix
- self.result[3][...] = (1.0 * self.field2[1][np.roll(ind,+2),...] -
- 8.0 * self.field2[1][np.roll(ind,+1),...] +\
- 8.0 * self.field2[1][np.roll(ind,-1),...] -\
- 1.0 * self.field2[1][np.roll(ind,-2),...]) / (12. * mesh_size[0])
-
- self.result[4][...] = (1.0 * self.field2[1][:,np.roll(ind,+2),:] -\
- 8.0 * self.field2[1][:,np.roll(ind,+1),:] +\
- 8.0 * self.field2[1][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[1][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
-
- self.result[5][...] = (1.0 * self.field2[1][...,np.roll(ind,+2)] -\
- 8.0 * self.field2[1][...,np.roll(ind,+1)] +\
- 8.0 * self.field2[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[1][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
-
-# Third line of Grad U matrix
- self.result[6][...] = (1.0 * self.field2[2][np.roll(ind,+2),...] -
- 8.0 * self.field2[2][np.roll(ind,+1),...] +\
- 8.0 * self.field2[2][np.roll(ind,-1),...] -\
- 1.0 * self.field2[2][np.roll(ind,-2),...]) / (12. * mesh_size[0])
-
- self.result[7][...] = (1.0 * self.field2[2][:,np.roll(ind,+2),:] -\
- 8.0 * self.field2[2][:,np.roll(ind,+1),:] +\
- 8.0 * self.field2[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field2[2][:,np.roll(ind,-2),:]) / (12. * mesh_size[1])
-
- self.result[8][...] = (1.0 * self.field2[2][...,np.roll(ind,+2)] -\
- 8.0 * self.field2[2][...,np.roll(ind,+1)] +\
- 8.0 * self.field2[2][...,np.roll(ind,-1)] -\
- 1.0 * self.field2[2][...,np.roll(ind,-2)]) / (12. * mesh_size[2])
+ maxgersh[0] = max(maxstr1,maxstr2,maxstr3)
+ maxgersh[1] = max(maxadv1,maxadv2,maxadv3)
+ result = np.concatenate((np.array([temp1, temp2, temp3]),np.array([temp4, temp5, temp6]),np.array([temp7, temp8, temp9])))
+ return result, maxgersh
elif self.choice == 'curl':
@@ -322,35 +207,39 @@ class DifferentialOperator_d(DiscreteOperator):
temp1 = (1.0 * self.field1[2][:,np.roll(ind,+2),:] -
8.0 * self.field1[2][:,np.roll(ind,+1),:] +\
8.0 * self.field1[2][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[2][:,np.roll(ind,-2),:] ) / (12. * mesh_size[1])
+ 1.0 * self.field1[2][:,np.roll(ind,-2),:] ) / (12. * self.meshSize[1])
temp2 = (1.0 * self.field1[1][...,np.roll(ind,+2)] -\
8.0 * self.field1[1][...,np.roll(ind,+1)] +\
8.0 * self.field1[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[1][...,np.roll(ind,-2)] ) / (12. * mesh_size[2])
- self.result[0]= temp1 - temp2
+ 1.0 * self.field1[1][...,np.roll(ind,-2)] ) / (12. * self.meshSize[2])
+ tmp1 = temp1 - temp2
+
+ temp1 = (1.0 * self.field1[0][...,np.roll(ind,+2)] -\
+ 8.0 * self.field1[0][...,np.roll(ind,+1)] +\
+ 8.0 * self.field1[0][...,np.roll(ind,-1)] -\
+ 1.0 * self.field1[0][...,np.roll(ind,-2)] ) / (12. * self.meshSize[2])
- temp1 = (1.0 * self.field1[2][np.roll(ind,+2),...] -
+ temp2 = (1.0 * self.field1[2][np.roll(ind,+2),...] -
8.0 * self.field1[2][np.roll(ind,+1),...] +\
8.0 * self.field1[2][np.roll(ind,-1),...] -\
- 1.0 * self.field1[2][np.roll(ind,-2),...] ) / (12. * mesh_size[0])
+ 1.0 * self.field1[2][np.roll(ind,-2),...] ) / (12. * self.meshSize[0])
- temp2 = (1.0 * self.field1[0][...,np.roll(ind,+2)] -\
- 8.0 * self.field1[0][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[0][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[0][...,np.roll(ind,-2)] ) / (12. * mesh_size[2])
- self.result[1]= temp1 - temp2
+ tmp2 = temp1 - temp2
- temp1 = (1.0 * self.field1[1][...,np.roll(ind,+2)] -
- 8.0 * self.field1[1][...,np.roll(ind,+1)] +\
- 8.0 * self.field1[1][...,np.roll(ind,-1)] -\
- 1.0 * self.field1[1][...,np.roll(ind,-2)] ) / (12. * mesh_size[0])
+ temp1 = (1.0 * self.field1[1][np.roll(ind,+2),...] -
+ 8.0 * self.field1[1][np.roll(ind,+1),...] +\
+ 8.0 * self.field1[1][np.roll(ind,-1),...] -\
+ 1.0 * self.field1[1][np.roll(ind,-2),...] ) / (12. * self.meshSize[0])
temp2 = (1.0 * self.field1[0][:,np.roll(ind,+2),:] -\
8.0 * self.field1[0][:,np.roll(ind,+1),:] +\
8.0 * self.field1[0][:,np.roll(ind,-1),:] -\
- 1.0 * self.field1[0][:,np.roll(ind,-2),:] ) / (12. * mesh_size[1])
- self.result[2]= temp1 - temp2
+ 1.0 * self.field1[0][:,np.roll(ind,-2),:] ) / (12. * self.meshSize[1])
+ tmp3 = temp1 - temp2
+ result = np.concatenate((np.array([tmp1]), np.array([tmp2]), np.array([tmp3])))
+ return result
+
else:
raise ValueError("Unknown differiential operator: " + str(self.choice))
diff --git a/HySoP/hysop/operator/fct2op.py b/HySoP/hysop/operator/fct2op.py
new file mode 100644
index 0000000000000000000000000000000000000000..485c57d9789250b0b9fe6cf0629dba33a1624976
--- /dev/null
+++ b/HySoP/hysop/operator/fct2op.py
@@ -0,0 +1,62 @@
+# -*- coding: utf-8 -*-
+"""
+@package operator
+Convert differential operator for stretching into a function to Discrete time integration
+"""
+from ..constants import *
+from .discrete import DiscreteOperator
+from .differentialOperator_d import DifferentialOperator_d
+import numpy as np
+
+class Fct2Op(object):
+ """
+ Fct2Op
+
+ make the link between differentialOperator_d and the integration time scheme (Euler,RK,...)
+ """
+
+ def __init__(self, field1, field2, choice=None, resolution = None , meshSize = None):
+ """
+ Constructor.
+ store field2 to work with it in the apply function
+
+ @param field1 : vorticity field
+ @param field2 : velocity field
+ @param choice : can be 'gradU' for w.Du or divWU for div(wu)
+ @param resolution : resolution of the integration domain (array of number of points)
+ @param meshSize : array of the space step
+ """
+ self.field1 = field1
+ self.field2 = field2
+ self.choice = choice
+ self.resolution = resolution
+ self.meshSize = meshSize
+
+ def apply(self, t , y):
+ """
+ Apply operator.
+
+ We used Fct2Op.apply like a def function
+ """
+ if self.choice == 'gradU' :
+ #Make Scalar product
+ temp0 = self.field2[0][:,:,:]* y[0][:,:,:] +self.field2[1][:,:,:]* y[1][:,:,:]+self.field2[2][:,:,:]* y[2][:,:,:]
+ temp1 = self.field2[3][:,:,:]* y[0][:,:,:] +self.field2[4][:,:,:]* y[1][:,:,:]+self.field2[5][:,:,:]* y[2][:,:,:]
+ temp2 = self.field2[6][:,:,:]* y[0][:,:,:] +self.field2[7][:,:,:]* y[1][:,:,:]+self.field2[8][:,:,:]* y[2][:,:,:]
+ result = np.concatenate((np.array([temp0]),np.array([temp1]),np.array([temp2])))
+ return result
+ if self.choice == 'divWU' :
+ result = DifferentialOperator_d(y, self.field2, self.choice, self.resolution , self.meshSize).apply()
+ return result
+ else :
+ print 'choice', choice , 'is not define'
+
+
+ def __str__(self):
+ s = "fct2Op(DiscreteOperator). " + DiscreteOperator.__str__(self)
+ return s
+
+if __name__ == "__main__":
+ print __doc__
+ print "- Provided class : Fct2Op"
+ print Fct2Op.__doc__
diff --git a/HySoP/hysop/operator/stretching.py b/HySoP/hysop/operator/stretching.py
index 21640cd26ed8b32aad4020e55123bd1664b58765..55d8fc1105b23663b8e157c61040224039937c74 100755
--- a/HySoP/hysop/operator/stretching.py
+++ b/HySoP/hysop/operator/stretching.py
@@ -6,6 +6,7 @@ Penalization operator representation
"""
from .continuous import ContinuousOperator
from .stretching_d import Stretching_d
+from ..particular_solvers.integrator.runge_kutta3 import RK3
class Stretching(ContinuousOperator):
@@ -28,18 +29,18 @@ class Stretching(ContinuousOperator):
self.vorticity = vorticity
self.addVariable([velocity, vorticity])
- def discretize(self, idVelocityD=0, idVorticityD=0, res_vorticity=None, method=None): ## rajouter un attribut méthode (avec RK4, RK2 etc pour integration en temps) ?
+ def discretize(self, idVelocityD=0, idVorticityD=0, propertyOp='divConservation', method=RK3):
"""
Stretching operator discretization method.
Create a discrete Stretching operator from given specifications.
@param idVelocityD : Index of velocity discretisation to use.
@param idVorticityD : Index of vorticity discretisation to use.
- @param res_vorticity : result stretch (vorticity).
- @param method : the method to use.
+ @param propertyOp : choice of the property garantied by the stretching Op ( divConservation or massConservation)
+ @param method : the method to use. (Euler, RK, ..) default : RK3
"""
if self.discreteOperator is None:
- self.discreteOperator = Stretching_d(self, idVelocityD, idVorticityD, res_vorticity, method)
+ self.discreteOperator = Stretching_d(self, idVelocityD, idVorticityD, propertyOp, method)
def __str__(self):
"""ToString method"""
diff --git a/HySoP/hysop/operator/stretching_d.py b/HySoP/hysop/operator/stretching_d.py
index 0271e9b8c2401beb111905620a6eab4267bb1a48..ada0a269cd833715331aa7d9c993011a40948afb 100755
--- a/HySoP/hysop/operator/stretching_d.py
+++ b/HySoP/hysop/operator/stretching_d.py
@@ -7,14 +7,16 @@ from ..constants import *
from .discrete import DiscreteOperator
from .differentialOperator import DifferentialOperator
from .differentialOperator_d import DifferentialOperator_d
+from .fct2op import Fct2Op
from ..particular_solvers.integrator.euler import Euler
-from ..particular_solvers.integrator.runge_kutta2stretching import RK2Stretch
-from ..particular_solvers.integrator.runge_kutta3stretching import RK3Stretch
-from ..particular_solvers.integrator.runge_kutta4stretching import RK4Stretch
+from ..particular_solvers.integrator.runge_kutta2 import RK2
+from ..particular_solvers.integrator.runge_kutta3 import RK3
+from ..particular_solvers.integrator.runge_kutta4 import RK4
from ..particular_solvers.integrator.integrator import ODESolver
import numpy as np
import numpy.testing as npt
import time
+import sys
class Stretching_d(DiscreteOperator):
@@ -23,7 +25,7 @@ class Stretching_d(DiscreteOperator):
DiscreteOperator.DiscreteOperator specialization.
"""
- def __init__(self, stretch, idVelocityD=0, idVorticityD=0, res_vorticity=None, method=None):
+ def __init__(self, stretch, idVelocityD=0, idVorticityD=0, propertyOp='divConservation', method=None):
"""
Constructor.
Create a Stretching operator on a given continuous domain.
@@ -37,16 +39,10 @@ class Stretching_d(DiscreteOperator):
self.velocity = stretch.velocity.discreteField[idVelocityD]
## Vorticity.
self.vorticity = stretch.vorticity.discreteField[idVorticityD]
- ## Result vorticity on the grid
- self.res_vorticity = res_vorticity.discreteField[idVorticityD]
- ## Stretching result on the grid
- self.stretchTermResult = [np.zeros((res_vorticity.discreteField[idVorticityD].topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(res_vorticity.discreteField[idVorticityD].topology.dim)]
- ## input fields
- self.input = [self.velocity, self.vorticity]
- ## output fields
- self.output = [self.res_vorticity]
+# ## input fields
+# self.input = [self.velocity, self.vorticity]
+ self.propertyOp = propertyOp
self.method = method
- self.maxgersh = np.zeros(1, dtype=PARMES_REAL, order=ORDER)
## self.name = "stretching"
def apply(self, t, dt):
@@ -57,8 +53,8 @@ class Stretching_d(DiscreteOperator):
@param dt : time step.
Stretching algorithm:
- @li 1. discretization of the divergence operator (div wu) : \n
- @li 2. Time integration. Performs a method choose to resolove the equation dw/dt = div(wu).
+ @li 1. discretization of the divergence operator (div wu) or w grad(u) : \n
+ @li 2. Time integration. Performs a method choose to resolove the equation dw/dt = div(wu) or dw/dt = w grad(u).
-define the function of ODEsolver as the result of the first step
-integrate with the chosen method (RK4. RK2. euler)
"""
@@ -66,39 +62,41 @@ class Stretching_d(DiscreteOperator):
if self.method is not None:
- ## Space discretisation of div (wu)
- self.stretchTerm = DifferentialOperator(self.stretch.vorticity, self.stretch.velocity, choice='divConservation')
- self.stretchTerm.discretize(self.vorticity, self.velocity, self.stretchTermResult, self.maxgersh)
- self.stretchTerm.discreteOperator.apply()
-
- # Time integration of dw/dt = div(wu)
-
- ndt = self.stabilityTest(dt, self.maxgersh, self.method)
- dtstretch = dt/float(ndt)
-
- for i in range (ndt) :
- if Euler == self.method :
- print "Euler"
- methodInt=self.method(self.stretchTerm.discreteOperator)
- for i in range(self.vorticity.topology.dim):
- self.res_vorticity[i][...] = methodInt.integrate(self.vorticity, self.stretchTermResult, t, dtstretch, direction=i)
-
- if RK2Stretch == self.method :
- print "RK2"
- methodInt=self.method(self.stretchTerm.discreteOperator)
- self.res_vorticity = methodInt.integrate(self.stretchTerm, self.stretchTermResult, self.vorticity, self.velocity, self.res_vorticity, t, dtstretch, self.vorticity.topology.dim)
-
- if RK3Stretch == self.method :
- print "RK3"
- methodInt=self.method(self.stretchTerm.discreteOperator)
- self.res_vorticity = methodInt.integrate(self.stretchTerm, self.stretchTermResult, self.vorticity, self.velocity, self.res_vorticity, t, dtstretch, self.vorticity.topology.dim)
-
- if RK4Stretch == self.method :
- print "RK4"
- methodInt=self.method(self.stretchTerm.discreteOperator)
- self.res_vorticity = methodInt.integrate(self.stretchTerm, self.stretchTermResult, self.vorticity, self.velocity, self.res_vorticity, t, dtstretch, self.vorticity.topology.dim)
-
- self.vorticity.data= np.copy(self.res_vorticity.data)
+
+ if self.propertyOp == 'divConservation' :
+
+ ## Space discretisation of grad(u)
+ self.stretchOp = DifferentialOperator_d(self.vorticity, self.velocity, choice='gradU', resolution=self.vorticity.topology.mesh.resolution, meshSize=self.vorticity.topology.mesh.size)
+ gradResult, maxgersh = self.stretchOp.apply()
+
+ ## Determination of Stretching time step (stability condition)
+ self.ndt = self.stabilityTest(dt, maxgersh[0], self.method)
+ self.dtstretch = dt/float(self.ndt)
+
+ print "Maxgersh", maxgersh
+ print "dtStretch, ndt", self.dtstretch, self.ndt
+ print "dtAdvec", 1./maxgersh[1]
+
+ ## fct2Op
+ self.fonction = Fct2Op(self.vorticity.data, gradResult,choice = 'gradU', resolution=self.vorticity.topology.mesh.resolution, meshSize=self.vorticity.topology.mesh.size)
+
+ elif self.propertyOp == 'massConservation' :
+ ## Determination of Stretching time step (stability condition)
+ self.ndt = 1
+ self.dtstretch = dt
+
+ ## fct2Op
+ self.fonction = Fct2Op(self.vorticity, self.velocity,choice = 'divWU', resolution=self.vorticity.topology.mesh.resolution, meshSize=self.vorticity.topology.mesh.size)
+
+ else:
+ raise ValueError("Unknown Stretching Operator property: " + str(self.propertyOp))
+
+ # Time integration
+ for i in range (self.ndt) :
+ methodInt=self.method(self.fonction)
+ self.res_vorticity = methodInt.integrate(self.fonction.apply, t, self.dtstretch, self.vorticity.data)
+
+ self.vorticity.data= np.copy(self.res_vorticity)
self.compute_time = time.time() - self.compute_time
self.total_time += self.compute_time
return self.compute_time
@@ -109,17 +107,17 @@ class Stretching_d(DiscreteOperator):
if method == Euler :
cststretch = 2.0
dtstretch = min(dtstretch, cststretch/maxi)
- if method == RK2Stretch :
+ if method == RK2 :
cststretch = 2.0
dtstretch = min(dtstretch, cststretch/maxi)
-# if method == RK3Stretch :
-# cststretch = 2.5127
-# dtstretch = min(dtstretch, cststretch/maxi)
- if method == RK4Stretch :
+ if method == RK3 :
+ cststretch = 2.5127
+ dtstretch = min(dtstretch, cststretch/maxi)
+ if method == RK4 :
cststretch = 2.7853
dtstretch = min(dtstretch, cststretch/maxi)
if dt == dtstretch :
- print "dt, ndt , dtstretch, upperbound", dt, int(dt/dtstretch), dt/float(int(dt/dtstretch)), 6./maxi
+ print "dt, ndt , dtstretch, upperbound", dt, int(dt/dtstretch), dt/float(int(dt/dtstretch)), cststretch/maxi
return int(dt/dtstretch)
else :
print "dt, ndt , dtstretch, upperbound", dt, int(dt/dtstretch)+1, dt/float(int(dt/dtstretch)+1), cststretch/maxi
diff --git a/HySoP/hysop/particular_solvers/basic.py b/HySoP/hysop/particular_solvers/basic.py
index c7300e269e4b49e0582ac469977cbb2500d9bba7..e694ecabc75f16c57bde3472abffccefed8053f3 100644
--- a/HySoP/hysop/particular_solvers/basic.py
+++ b/HySoP/hysop/particular_solvers/basic.py
@@ -7,9 +7,8 @@ from solver import Solver
#from integrator.runge_kutta4 import RK4
from .integrator.integrator import ODESolver
from .integrator.euler import Euler
-from .integrator.runge_kutta2stretching import RK2Stretch
-from .integrator.runge_kutta4stretching import RK4Stretch
from .integrator.runge_kutta2 import RK2
+from .integrator.runge_kutta3 import RK3
from .integrator.runge_kutta4 import RK4
from ..fields.continuous import ContinuousField
from ..operator.transport import Transport
@@ -52,8 +51,7 @@ class ParticleSolver(Solver):
## Advection operator of the problem. Advection need special treatment in particle methods.
self.advection = None
## ODE Solver method.
-# self.ODESolver = ODESolver
- self.ODESolver = Euler
+ self.ODESolver = ODESolver
## Interpolation method.
self.InterpolationMethod = InterpolationMethod
## Remeshing method.
@@ -86,9 +84,8 @@ class ParticleSolver(Solver):
Initialize a particle method solver regarding the problem.
"""
self.p_position = ContinuousField(domain=self.problem.domains[0], name='ParticlePosition')
- self.p_vorticity = ContinuousField(domain=self.problem.domains[0], name='VortiStretch',vector=True)
self.p_scalar = ContinuousField(domain=self.problem.domains[0], name='ParticleScalar')
- self.problem.addVariable([self.p_position, self.p_scalar, self.p_vorticity])
+ self.problem.addVariable([self.p_position, self.p_scalar])
## Discretise domains
for d in self.problem.domains:
d.discretize(self.problem.topology.resolution)
@@ -104,7 +101,7 @@ class ParticleSolver(Solver):
if op is self.advection:
op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
elif op is self.stretch:
- op.discretize(res_vorticity=self.p_vorticity, method=self.ODESolver)
+ op.discretize(method=self.ODESolver)
## Velocity is only program on gpu for instance
#elif op is self.velocity:
# op.discretize(result_position=self.p_position, result_scalar=self.p_scalar, method=self.ODESolver)
diff --git a/HySoP/hysop/particular_solvers/integrator/euler.py b/HySoP/hysop/particular_solvers/integrator/euler.py
index 7fbca88149e49f45a258ac21189a4c3091e57d1e..0a656e695ee4e5e257d1e3d9912493d3f748736d 100755
--- a/HySoP/hysop/particular_solvers/integrator/euler.py
+++ b/HySoP/hysop/particular_solvers/integrator/euler.py
@@ -28,19 +28,7 @@ class Euler(ODESolver):
"""
ODESolver.__init__(self, f)
-# def integrate(self, y, f, t, dt, direction):
-# """
-# Integration step for forward Euler method.
-# y(t_{n+1}) = y(t_n) + dt*f(y(t_n))
-
-# @param y : position at time t.
-# @param t : current time.
-# @param dt : time step.
-# @return y at t+dt.
-# """
-# return y[direction][...] +dt * f[direction][...]
-
- def integrate(self, f, t , dt , u, topo ):
+ def integrate(self, f, t, dt, y):
"""
Integration step for Euler Method.
up = f[t, u(t,x,y,z)].
@@ -50,7 +38,7 @@ class Euler(ODESolver):
@param t : current time.
@param dt : time step.
# """
- result = u[:,:,:,:] + dt * np.asarray(f(t,u[:,:,:,:]))
+ result = y[:,:,:,:] + dt * np.asarray(f(t,y[:,:,:,:]))
return result
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
index 98982cf82b2afcbd7630ed0d3815ca8592ff9a10..6bff109502efb91c1a665bc4be1e7a326c1fc8fa 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta2.py
@@ -21,38 +21,24 @@ class RK2(ODESolver):
@param f function f(t,y) : Right hand side of the equation to solve.
@param dim : dimensions
- @param conditions : function to apply boundary conditions.
"""
ODESolver.__init__(self, f)
-# def integrate(self, y, fy, yp, t, dt, dir):
-# """
-# Integration step for RK2 Method.
-# yp= y + dt*y'[y+dt/2*y'(y)].
-
-# @param y : position at time t.
-# @param fy : y'(y) value at time t.
-# @param yp : result position at time t + dt
-# @param t : current time.
-# @param dt : time step.
-# """
-# pass
-
- def integrate(self, f, t , dt , u, topo ):
+ def integrate(self, f, t , dt , y):
"""
Integration step for RK2 Method.
- up = f[t, u(t,x,y,z)].
+ up = f[t, Y(t,x,y,z)].
- @param u : function of position at time t.
+ @param y : function of position at time t.
@param f : function of function of position at time t.
@param t : current time.
@param dt : time step.
"""
-#K1 = dt *f[t,u(t,x,y,z)]
- K1 = dt * np.asarray(f(t,u))
-#result = u(t,x,y,z)] + dt *f[t + dt/2 , u(t,x,y,z) + k1/2.]
- result = u + dt *np.asarray(f(t+dt/2.,u + K1/2.))
+#K1 = dt *f[t,y(t,x,y,z)]
+ K1 = dt * np.asarray(f(t,y))
+#result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+ result = y + dt *np.asarray(f(t+dt/2.,y + K1/2.))
return result
def __str__(self):
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
index 25e7b282c9a30d6a0ce4d45eaa3a4f7f43bd814a..bf2c868071dc7cee795783696720c32c48334b43 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta2stretching.py
@@ -45,14 +45,13 @@ class RK2Stretch(ODESolver):
@param t : current time.
@param dt : time step.
"""
- print "ehoh bis"
maxgersh = np.zeros(1, dtype=PARMES_REAL, order=ORDER)
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=K1, field2=field2, result=resultTmp, choice='divConservation', maxgersh=maxgersh)
+ f.discreteOperator = DifferentialOperator_d(field1=K1, field2=field2, result=resultTmp, choice='divConservation', maxgersh=maxgersh)
f.discreteOperator.apply()
res_vorticity[:,:,:] = field1[:,:,:] + resultTmp[:,:,:,:] * dt
return res_vorticity
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
index 6438426eff0e61422b1395ca29b9bc1877892fde..dc1d95beab3c67826455acd1df6eae5348c12224 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta3.py
@@ -25,37 +25,25 @@ class RK3(ODESolver):
"""
ODESolver.__init__(self, f)
-# def integrate(self, y, fy, yp, t, dt, dir):
-# """
-# Integration step for RK3 Method.
-# yp= y + dt*y'[y+dt/2*y'(y)].
-# @param y : position at time t.
-# @param fy : y'(y) value at time t.
-# @param yp : result position at time t + dt
-# @param t : current time.
-# @param dt : time step.
-# """
-# pass
-
- def integrate(self, f, t , dt , u, topo ):
+ def integrate(self, f, t , dt , y):
"""
Integration step for RK3 Method.
- up = f[t, u(t,x,y,z)].
+ up = f[t, Y(t,x,y,z)].
- @param u : function of position at time t.
+ @param y : function of position at time t.
@param f : function of function of position at time t.
@param t : current time.
@param dt : time step.
"""
-#K1 = dt *f[t,u(t,x,y,z)]
- K1 = dt * np.asarray(f(t,u))
-#result = u(t,x,y,z)] + dt *f[t + dt/2 , u(t,x,y,z) + k1/2.]
- K2 = dt *np.asarray(f(t+dt/3.,u + K1[:,:,:,:]/3.))
- K3 = dt *np.asarray(f(t+dt*2./3.,u + K2[:,:,:,:]*2./3.))
-#result = u(t,x,y,z)] + 1/4 [ k1 + 3 k3]
- result = u + 1./4.*(K1+3.*K3)
+#K1 = dt *f[t,y(t,x,y,z)]
+ K1 = dt * np.asarray(f(t,y))
+#result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+ K2 = dt *np.asarray(f(t+dt/3.,y + K1[:,:,:,:]/3.))
+ K3 = dt *np.asarray(f(t+dt*2./3.,y + K2[:,:,:,:]*2./3.))
+#result = y(t,x,y,z)] + 1/4 [ k1 + 3 k3]
+ result = y + 1./4.*(K1+3.*K3)
return result
def __str__(self):
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta3stretching.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta3stretching.py
index 789659b3c2a8e4507334033a398300d0d70da3d8..9d9c62ecff9912aa3d81384ff091f343c6cfb5e6 100644
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta3stretching.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta3stretching.py
@@ -47,10 +47,10 @@ class RK3Stretch(ODESolver):
K1=np.asarray(fd)
K2=np.asarray([np.zeros((field2.topology.resolution), dtype=PARMES_REAL, order=ORDER) for d in xrange(field2.topology.dim)])
K1[:,:,:,:] = field1[:,:,:] + dt * K1[:,:,:,:]
- f.discreteOperator = DifferentialOperator_d(f, K1, field2, K2, choice='divConservation',maxgersh=maxgersh)
+ f.discreteOperator = DifferentialOperator_d(K1, field2, K2, choice='divConservation',maxgersh=maxgersh)
f.discreteOperator.apply()
K2[:,:,:,:] = 3./4. * field1[:,:,:] + 1./4. * (K1[:,:,:,:] + dt * K2[:,:,:,:])
- f.discreteOperator = DifferentialOperator_d(f, K2, field2, res_vorticity, choice='divConservation',maxgersh=maxgersh)
+ f.discreteOperator = DifferentialOperator_d(K2, field2, res_vorticity, choice='divConservation',maxgersh=maxgersh)
f.discreteOperator.apply()
res_vorticity[:,:,:] = 1./3. * field1[:,:,:] + 2./3. * (K2[:,:,:,:] + dt * res_vorticity[:,:,:])
return res_vorticity
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
index 1bbedadc1b25c4eceecc7d6076076cdfb7749453..26bdba75037540614967e3da5b7ba0085143d90f 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta4.py
@@ -25,37 +25,24 @@ class RK4(ODESolver):
"""
ODESolver.__init__(self, f)
-# def integrate(self, y, fy, yp, t, dt, dir):
-# """
-# Integration step for RK4 Method.
-# yp= y + dt*y'[y+dt/2*y'(y)].
-
-# @param y : position at time t.
-# @param fy : y'(y) value at time t.
-# @param yp : result position at time t + dt
-# @param t : current time.
-# @param dt : time step.
-# """
-# pass
-
- def integrate(self, f, t , dt , u, topo ):
+ def integrate(self, f, t , dt , y):
"""
Integration step for RK4 Method.
- up = f[t, u(t,x,y,z)].
+ up = f[t, Y(t,x,y,z)].
- @param u : function of position at time t.
+ @param y : function of position at time t.
@param f : function of function of position at time t.
@param t : current time.
@param dt : time step.
"""
-#K1 = dt *f[t,u(t,x,y,z)]
- K1 = dt * np.asarray(f(t,u))
- K2 = dt *np.asarray(f(t+dt/2.,u + K1[:,:,:,:]/2.))
- K3 = dt *np.asarray(f(t+dt/2.,u + K2[:,:,:,:]/2.))
- K4 = dt *np.asarray(f(t+dt,u + K3[:,:,:,:]))
-#result = u(t,x,y,z)] + dt *f[t + dt/2 , u(t,x,y,z) + k1/2.]
- result = u + 1./6. * (K1[:,:,:,:] + 2.*K2[:,:,:,:] + 2.*K3[:,:,:,:]+ K4[:,:,:,:])
+#K1 = dt *f[t,y(t,x,y,z)]
+ K1 = dt * np.asarray(f(t,y))
+ K2 = dt *np.asarray(f(t+dt/2.,y + K1[:,:,:,:]/2.))
+ K3 = dt *np.asarray(f(t+dt/2.,y + K2[:,:,:,:]/2.))
+ K4 = dt *np.asarray(f(t+dt,y + K3[:,:,:,:]))
+#result = y(t,x,y,z)] + dt *f[t + dt/2 , y(t,x,y,z) + k1/2.]
+ result = y + 1./6. * (K1[:,:,:,:] + 2.*K2[:,:,:,:] + 2.*K3[:,:,:,:]+ K4[:,:,:,:])
return result
def __str__(self):
diff --git a/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py b/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
index b16238fbddedc59cda83b09ba80623e9c4e550aa..6427adefaff4a613feab97de4ae38494ae1f901a 100755
--- a/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
+++ b/HySoP/hysop/particular_solvers/integrator/runge_kutta4stretching.py
@@ -51,15 +51,15 @@ class RK4Stretch(ODESolver):
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='divConservation',maxgersh=maxgersh)
+ f.discreteOperator = DifferentialOperator_d(resultTmp , field2, K2, choice='divConservation',maxgersh=maxgersh)
f.discreteOperator.apply()
K2[:,:,:,:] = K2[:,:,:,:] * dt
resultTmp[:,:,:,:] = K2[:,:,:,:] * 0.5 + field1[:,:,:]
- f.discreteOperator = DifferentialOperator_d(f, resultTmp, field2, K3, choice='divConservation',maxgersh=maxgersh)
+ f.discreteOperator = DifferentialOperator_d(resultTmp, field2, K3, choice='divConservation',maxgersh=maxgersh)
f.discreteOperator.apply()
K3[:,:,:,:] = K3[:,:,:,:] * dt
resultTmp[:,:,:,:] = K3[:,:,:,:] + field1[:,:,:]
- f.discreteOperator = DifferentialOperator_d(f, resultTmp, field2, K4, choice='divConservation',maxgersh=maxgersh)
+ f.discreteOperator = DifferentialOperator_d(resultTmp, field2, K4, choice='divConservation',maxgersh=maxgersh)
f.discreteOperator.apply()
K4[:,:,:,:] = K4[:,:,:,:] * dt
res_vorticity[:,:,:] = field1[:,:,:] + 1./6. * (K1[:,:,:,:] + 2.0 * K2[:,:,:,:] + 2.0 * K3[:,:,:,:] + K4[:,:,:,:])
diff --git a/HySoP/hysop/problem/problem.py b/HySoP/hysop/problem/problem.py
index bceb93e6995eb14c8901f7f63af80685647dfc54..fccd368596785707ef19b95859e0f48154a0059c 100644
--- a/HySoP/hysop/problem/problem.py
+++ b/HySoP/hysop/problem/problem.py
@@ -93,7 +93,7 @@ class Problem():
self.initSolver()
self.io.step()
while not self.timer.end:
- print "==== Iteration : {0:3d} t={1:6.2f} ====".format(self.timer.ite + 1, self.timer.t + self.timer.dt)
+ print "==== Iteration : {0:3d} t={1:6.3f} ====".format(self.timer.ite + 1, self.timer.t + self.timer.dt)
for op in self.operators:
op.apply(self.timer.t, self.timer.dt)
self.timer.step()
diff --git a/HySoP/hysop/test/data/Fields_sav0_U.data b/HySoP/hysop/test/data/Fields_sav0_U.data
new file mode 100644
index 0000000000000000000000000000000000000000..4406cbd096d1b19e3e1f467467f3821b6a945209
Binary files /dev/null and b/HySoP/hysop/test/data/Fields_sav0_U.data differ
diff --git a/HySoP/hysop/test/data/Fields_sav0_V.data b/HySoP/hysop/test/data/Fields_sav0_V.data
new file mode 100644
index 0000000000000000000000000000000000000000..3beca511fc042cf0897a96300f94661ad834494a
Binary files /dev/null and b/HySoP/hysop/test/data/Fields_sav0_V.data differ
diff --git a/HySoP/hysop/test/data/Fields_sav0_W.data b/HySoP/hysop/test/data/Fields_sav0_W.data
new file mode 100644
index 0000000000000000000000000000000000000000..d2637b164b546619c9f66a98695c3900d30cfee7
Binary files /dev/null and b/HySoP/hysop/test/data/Fields_sav0_W.data differ
diff --git a/HySoP/hysop/test/test_operator/test_CondStability.py b/HySoP/hysop/test/test_operator/test_CondStability.py
new file mode 100755
index 0000000000000000000000000000000000000000..652751e6f4fb509b6481db14f0e69ff2be0c1f2a
--- /dev/null
+++ b/HySoP/hysop/test/test_operator/test_CondStability.py
@@ -0,0 +1,141 @@
+# -*- coding: utf-8 -*-
+import time
+from parmepy.operator.differentialOperator_d import DifferentialOperator_d
+from parmepy.particular_solvers.integrator.euler import Euler
+from parmepy.particular_solvers.integrator.runge_kutta2 import RK2
+from parmepy.particular_solvers.integrator.runge_kutta3 import RK3
+from parmepy.particular_solvers.integrator.runge_kutta4 import RK4
+import parmepy as pp
+from parmepy.constants import *
+import numpy as np
+from math import *
+import unittest
+import sys
+import struct
+import array
+
+
+class test_CondStability(unittest.TestCase):
+ """
+ Condition Stability test class
+ """
+
+ def vitesse(self,x, y, z):
+ vx = 1.
+ vy = 1.
+ vz = 1.
+ return vx, vy, vz
+
+ def vorticite(self,x, y, z):
+ wx = 1.
+ wy = 1.
+ wz = 1.
+ return wx, wy, wz
+
+ def scalaire(self,x, y, z):
+ if x < 0.5 and y < 0.5 and z < 0.5:
+ return 1.
+ else:
+ return 0.
+
+
+ def testCondStab(self):
+ # Parameters
+ nb = 128
+ timeStep = 0.09
+ finalTime = 0.09
+ self.t = 0.
+ t0 = time.time()
+
+ ## Domain
+ box = pp.Box(3, length=[1., 1., 1.], origin=[0., 0., 0.])
+
+ ## Fields
+ velo = pp.AnalyticalField(domain=box, formula=self.vitesse, name='Velocity', vector=True)
+ vorti = pp.AnalyticalField(domain=box, formula=self.vorticite, name='Vorticity', vector=True)
+
+ inputJBField = [np.zeros((128,128,128), dtype=PARMES_REAL, order=ORDER) for d in xrange(3)]
+
+ f1 = open('./parmepy/test/data/Fields_sav0_U.data','rb')
+ f2 = open('./parmepy/test/data/Fields_sav0_V.data','rb')
+ f3 = open('./parmepy/test/data/Fields_sav0_W.data','rb')
+
+ nbx = np.asarray(struct.unpack("i",f1.read(4)))
+ nby = np.asarray(struct.unpack("i",f1.read(4)))
+ nbz = np.asarray(struct.unpack("i",f1.read(4)))
+
+ binvalues = array.array('d')
+ binvalues.read(f1, nbx*nby*nbz)
+
+ data = np.array(binvalues, dtype=PARMES_REAL)
+ inputJBField[0] = np.reshape(data, (nbx,nby,nbz))
+ f1.close()
+
+ nbx = np.asarray(struct.unpack("i",f2.read(4)))
+ nby = np.asarray(struct.unpack("i",f2.read(4)))
+ nbz = np.asarray(struct.unpack("i",f2.read(4)))
+
+ binvalues = array.array('d')
+ binvalues.read(f2, nbx*nby*nbz)
+
+ data = np.array(binvalues, dtype=PARMES_REAL)
+ inputJBField[1] = np.reshape(data, (nbx,nby,nbz))
+ f2.close()
+
+ nbx = np.asarray(struct.unpack("i",f3.read(4)))
+ nby = np.asarray(struct.unpack("i",f3.read(4)))
+ nbz = np.asarray(struct.unpack("i",f3.read(4)))
+
+ binvalues = array.array('d')
+ binvalues.read(f3, nbx*nby*nbz)
+
+ data = np.array(binvalues, dtype=PARMES_REAL)
+ inputJBField[2] = np.reshape(data, (nbx,nby,nbz))
+ f3.close()
+
+ ## Operators
+ stretch = pp.Stretching(velo,vorti)
+
+ ## Solver creation (discretisation of objects is done in solver initialisation)
+ topo3D = pp.CartesianTopology(domain=box, resolution=[nbx[0], nby[0], nbz[0]], dim=3)
+
+ ##Problem
+ pb = pp.Problem(topo3D, [stretch])
+
+ ## Setting solver to Problem
+ pb.setSolver(finalTime, timeStep, solver_type='basic')
+ pb.solver.ODESolver= Euler#RK4# RK3# RK2#
+ pb.initSolver()
+
+
+ self.result = [np.ones((nbx, nby, nbz), dtype=PARMES_REAL, order=ORDER) for d in xrange(3)]
+ vortidata= [np.zeros((128,128,128), dtype=PARMES_REAL, order=ORDER) for d in xrange(3)]
+# print 'stretch', dir(stretch.variables), stretch.velocity.discreteField[0].data
+ stretch.velocity.discreteField[0].data = inputJBField
+
+ self.curl = DifferentialOperator_d(stretch.velocity.discreteField[0].data, vortidata, choice='curl', resolution=topo3D.mesh.resolution, meshSize=np.array([1./nbx[0], 1./nby[0], 1./nbz[0]]))
+# self.curl.discretize(velo, vorti, resolution=topo3D.mesh.resolution, meshSize=np.array(1./nbx[0], 1./nby[0], 1./nbz[0]))
+ stretch.vorticity.discreteField[0].data = self.curl.discreteOperator.apply()
+ t1 = time.time()
+
+ ## Solve problem
+ pb.solve()
+
+ tf = time.time()
+
+ print "\n"
+ print "Total time : ", tf - t0, "sec (CPU)"
+ print "Init time : ", t1 - t0, "sec (CPU)"
+ print "Solving time : ", tf - t1, "sec (CPU)"
+
+
+ def runTest(self):
+ self.testCondStab()
+
+def suite():
+ suite = unittest.TestSuite()
+ suite.addTest(unittest.makeSuite(test_CondStability))
+ return suite
+
+if __name__ == "__main__":
+ unittest.TextTestRunner(verbosity=2).run(suite())
diff --git a/HySoP/hysop/test/test_operator/test_Curl.py b/HySoP/hysop/test/test_operator/test_Curl.py
new file mode 100755
index 0000000000000000000000000000000000000000..c7131e1e87f8d33cec1fc60cc8920f4c5e501147
--- /dev/null
+++ b/HySoP/hysop/test/test_operator/test_Curl.py
@@ -0,0 +1,109 @@
+# -*- coding: utf-8 -*-
+import unittest
+
+#import parmepy
+
+from parmepy.operator.transport import *
+from parmepy.operator.continuous import *
+from parmepy.operator.differentialOperator import *
+from parmepy.operator.stretching import *
+from parmepy.fields.discrete import *
+from parmepy.fields.continuous import *
+from parmepy.fields.analytical import *
+from parmepy.domain.topology import *
+from parmepy.domain.box import *
+from parmepy.constants import *
+from parmepy.particular_solvers.basic import *
+from parmepy.particular_solvers.integrator.euler import *
+from parmepy.particular_solvers.solver import *
+import numpy as np
+import numpy.testing as npt
+import math
+
+
+class test_Curl(unittest.TestCase):
+ """
+ DiscreteVariable test class
+ """
+ def setUp(self):
+ self.e = 0.002 # Accepted error between result and analytical result
+ self.dim = 3
+ self.boxLength = [2.*np.pi, 2.*np.pi, 2.*np.pi]
+ self.boxMin = [ 0., 0., 0.]
+ self.nbPts = [32, 32, 32]
+ self.t = 0.
+ self.timeStep = 0.02
+ self.box = Box(dimension=self.dim,
+ length=self.boxLength,
+ origin=self.boxMin)
+
+ def testOperatorCurl(self):
+ # Continuous fields and operator declaration
+ self.velo = AnalyticalField(domain=self.box, formula=self.vitesse, name='Velocity', vector=True)
+ self.vorti = AnalyticalField(domain=self.box, formula=self.vorticite, name='Vorticity', vector=True)
+ # chercher cas test de tel sorte que le stretching serait periodique et analytique
+ self.curl = DifferentialOperator(self.velo, self.vorti, choice='curl')
+ # Topology definition / Fields and operator discretization
+ self.result = [np.ones((self.nbPts), dtype=PARMES_REAL, order=ORDER) for d in xrange(self.dim)]
+ self.topo3D = CartesianTopology(domain=self.box, resolution=self.nbPts, dim=self.dim)
+ self.vorti.discretize(self.topo3D)
+ self.velo.discretize(self.topo3D)
+ self.vorti.initialize()
+ self.velo.initialize()
+ self.curl.discretize(self.velo.discreteField[self.velo._fieldId], self.vorti.discreteField[self.vorti._fieldId], result=self.result)
+ self.curl.discreteOperator.apply()
+ self.FinalTime = 0.02
+ self.anal=np.vectorize(self.vorticite)(self.topo3D.mesh.coords[0], \
+ self.topo3D.mesh.coords[1], \
+ self.topo3D.mesh.coords[2])
+ # Comparison with analytical solution
+# print "max:",np.max(abs(self.anal[0]-self.result[0]))
+ npt.assert_array_less(abs(self.anal[0]-self.result[0]), self.e)
+# print "max:",np.max(abs(self.anal[1]-self.result[1]))
+ npt.assert_array_less(abs(self.anal[1]-self.result[1]), self.e)
+# print "max:",np.max(abs(self.anal[2]-self.result[2]))
+ npt.assert_array_less(abs(self.anal[2]-self.result[2]), self.e)
+
+ def vitesse(self, x, y, z):
+# amodul = np.cos(np.pi*self.t/3)
+# pix = np.pi*x
+# piy = np.pi*y
+# piz = np.pi*z
+# pi2x = 2.*pix
+# pi2y = 2.*piy
+# pi2z = 2.*piz
+# vx = 2.*np.sin(pix)*np.sin(pix)*np.sin(pi2y)*np.sin(pi2z)*amodul
+# vy = -np.sin(pi2x)*np.sin(piy)*np.sin(piy)*np.sin(pi2z)*amodul
+# vz = -np.sin(pi2x)*np.sin(piz)*np.sin(piz)*np.sin(pi2y)*amodul
+ vx = np.cos(y)
+ vy = np.cos(z)
+ vz = np.cos(x)
+ return vx, vy, vz
+
+ def vorticite(self, x, y, z):
+# amodul = np.cos(np.pi*self.t/3)
+# pix = np.pi*x
+# piy = np.pi*y
+# piz = np.pi*z
+# pi2x = 2.*pix
+# pi2y = 2.*piy
+# pi2z = 2.*piz
+# wx = 2.* np.pi * np.sin(pi2x) * amodul*( - np.cos(pi2y)*np.sin(piz)*np.sin(piz)+ np.sin(piy)*np.sin(piy)*np.cos(pi2z) )
+# wy = 2.* np.pi * np.sin(pi2y) * amodul*( 2.*np.cos(pi2z)*np.sin(pix)*np.sin(pix)+ np.sin(piz)*np.sin(piz)*np.cos(pi2x) )
+# wz = -2.* np.pi * np.sin(pi2z) * amodul*( np.cos(pi2x)*np.sin(piy)*np.sin(piy)+ np.sin(pix)*np.sin(pix)*np.cos(pi2y) )
+ wx = np.sin(z)
+ wy = np.sin(x)
+ wz = np.sin(y)
+ return wx, wy, wz
+
+ def runTest(self):
+ self.setUp()
+ self.testOperatorCurl()
+
+def suite():
+ suite = unittest.TestSuite()
+ suite.addTest(unittest.makeSuite(test_Curl))
+ return suite
+
+if __name__ == "__main__":
+ unittest.TextTestRunner(verbosity=2).run(suite())
diff --git a/HySoP/hysop/test/test_operator/test_Stretching.py b/HySoP/hysop/test/test_operator/test_Stretching.py
index 6d8ee863fa8b9f386a195da77bb899c3858c1d9e..b1588d92b54c9f8d5a903ff1641c595608f9d27f 100755
--- a/HySoP/hysop/test/test_operator/test_Stretching.py
+++ b/HySoP/hysop/test/test_operator/test_Stretching.py
@@ -1,9 +1,9 @@
# -*- coding: utf-8 -*-
import time
from parmepy.particular_solvers.integrator.euler import Euler
-from parmepy.particular_solvers.integrator.runge_kutta2stretching import RK2Stretch
-from parmepy.particular_solvers.integrator.runge_kutta3stretching import RK3Stretch
-from parmepy.particular_solvers.integrator.runge_kutta4stretching import RK4Stretch
+from parmepy.particular_solvers.integrator.runge_kutta2 import RK2
+from parmepy.particular_solvers.integrator.runge_kutta3 import RK3
+from parmepy.particular_solvers.integrator.runge_kutta4 import RK4
import parmepy as pp
import numpy as np
from math import *
@@ -81,7 +81,7 @@ class test_Stretching(unittest.TestCase):
## Setting solver to Problem
pb.setSolver(finalTime, timeStep, solver_type='basic')
- pb.solver.ODESolver= RK3Stretch# RK2Stretch# Euler#RK4Stretch#
+ pb.solver.ODESolver= RK4# RK3# RK2# Euler#
pb.initSolver()
t1 = time.time()