add liebermanIteration function

meteoinfo-lab/milconfig.xml

@@ -19,16 +19,16 @@
   </Path>
   <File>
     <OpenedFiles>
-      <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_3.py"/>
       <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_2.py"/>
       <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_1.py"/>
       <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_2.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_3.py"/>
     </OpenedFiles>
     <RecentFiles>
-      <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_3.py"/>
       <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_2.py"/>
       <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_1.py"/>
       <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_2.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_3.py"/>
     </RecentFiles>
   </File>
   <Font>
@@ -36,5 +36,5 @@
   </Font>
   <LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
   <Figure DoubleBuffering="true"/>
-  <Startup MainFormLocation="-7,0" MainFormSize="1398,806"/>
+  <Startup MainFormLocation="-7,0" MainFormSize="1493,844"/>
 </MeteoInfo>

meteoinfo-lab/pylib/mipylib/meteolib/calc/kinematics.py

@@ -55,20 +55,26 @@ def _grad_atmos(longitude, latitude, H):
     grad_y = grad_y / dy
     return grad_x, grad_y
 
-def velocity_potential(longitude, latitude, u, v, loop_max=1e10, epsilon=1e-10):
+def velocity_potential(longitude, latitude, u, v, loop_max=1e6, epsilon=1e-7, sor_index=0.2):
     """
     Calculate velocity potential using Richardson iterative method.
 
     Parameters
     ----------
     lon : (M, N) `array`
-        longitude array
+        longitude array.
     lat : (M, N) `array`
-        latitude array
+        latitude array.
     u : (M, N) `array`
-        x component of the wind
+        x component of the wind.
     v : (M, N) `array`
-        y component of the wind
+        y component of the wind.
+    loop_max : `int`
+        Maximum iteration loop number. Default is `1e6`.
+    epsilon : `float`
+        Minimum error value. Default is `1e-7`
+    sor_index : `float`
+        Super relaxation coefficient [0.2 - 0.5]. Default is `0.2`.
 
     Returns
     -------
@@ -81,8 +87,8 @@ def velocity_potential(longitude, latitude, u, v, loop_max=1e10, epsilon=1e-10):
     divh = _divh_atmos(longitude, latitude, u, v)  # vertical divergence
     dx2 = _dx_atmos(longitude, latitude)**2        # square of latitude gradient
     dy2 = _dy_atmos(latitude)**2                   # square of longitude gradient
-    phi = np.NDArray(MeteoMath.richardsonIteration(divh._array, dx2._array, dy2._array,
+    phi = np.NDArray(MeteoMath.liebermanIteration(divh._array, dx2._array, dy2._array,
-                                                   loop_max, epsilon))
+                                                   loop_max, epsilon, sor_index))
     phi = -phi
 
     # divergence wind
@@ -98,7 +104,7 @@ def velocity_potential(longitude, latitude, u, v, loop_max=1e10, epsilon=1e-10):
 
     return phi, Uphi, Vphi
 
-def stream_function(longitude, latitude, u, v, loop_max=int(1e10), epsilon=1e-10):
+def stream_function(longitude, latitude, u, v, loop_max=int(1e10), epsilon=1e-10, sor_index=0.2):
     """
     Calculate stream function using Richardson iterative method.
 
@@ -112,6 +118,12 @@ def stream_function(longitude, latitude, u, v, loop_max=int(1e10), epsilon=1e-10
         x component of the wind
     v : (M, N) `array`
         y component of the wind
+    loop_max : `int`
+        Maximum iteration loop number. Default is `1e6`.
+    epsilon : `float`
+        Minimum error value. Default is `1e-7`
+    sor_index : `float`
+        Super relaxation coefficient [0.2 - 0.5]. Default is `0.2`.
 
     Returns
     -------
@@ -124,8 +136,8 @@ def stream_function(longitude, latitude, u, v, loop_max=int(1e10), epsilon=1e-10
     curlz = _curlz_atmos(longitude, latitude, u, v)  # vorticity
     dx2 = _dx_atmos(longitude, latitude)**2        # square of latitude gradient
     dy2 = _dy_atmos(latitude)**2                   # square of longitude gradient
-    psi = np.NDArray(MeteoMath.richardsonIteration(curlz._array, dx2._array, dy2._array,
+    psi = np.NDArray(MeteoMath.liebermanIteration(curlz._array, dx2._array, dy2._array,
-                                                   loop_max, epsilon))
+                                                   loop_max, epsilon, sor_index))
     psi = -psi
 
     # vorticity wind

meteoinfo-math/src/main/java/org/meteoinfo/math/meteo/MeteoMath.java

@@ -1154,9 +1154,9 @@ public class MeteoMath {
                     // calculate residual
                     res.setDouble(i * n + j, (phi.getDouble((i + 1) * n + j) +
                             phi.getDouble((i - 1) * n + j) - 2 * phi.getDouble(i * n + j)) /
-                            dx2.getDouble(i * n + j) + (phi.getDouble(i * n + j + 1) +
+                            dy2.getDouble(i * n + j) + (phi.getDouble(i * n + j + 1) +
                             phi.getDouble(i * n + j - 1) - 2 * phi.getDouble(i * n + j)) /
-                            dy2.getDouble(i * n + j) + array.getDouble(i * n + j));
+                            dx2.getDouble(i * n + j) + array.getDouble(i * n + j));
                     // iterator
                     phi.setDouble(i * n + j, phi.getDouble(i * n + j) + res.getDouble(i * n + j) /
                             (2 / dx2.getDouble(i * n + j) + 2 / dy2.getDouble(i * n + j)));
@@ -1171,4 +1171,45 @@ public class MeteoMath {
         return phi;
     }
 
+    /**
+     * Lieberman's iteration method for calculating Poisson equation
+     *
+     * @param array Input array - 2D
+     * @param dx2 X gradient square
+     * @param dy2 Y gradient square
+     * @param loopMax Maximum loop number
+     * @param epsilon Minimum error value
+     * @param sorIndex Super relaxation coefficient
+     * @return Poisson equation output
+     */
+    public static Array liebermanIteration(Array array, Array dx2, Array dy2, long loopMax, double epsilon,
+                                           double sorIndex) {
+        int[] shape = array.getShape();
+        int m = shape[0];
+        int n = shape[1];
+        Array phi = ArrayUtil.zeros(shape, DataType.DOUBLE);
+        Array res = ArrayUtil.full(shape, -9999., DataType.DOUBLE);
+        for (int k=0; k < loopMax; k++) {
+            for (int i = 1; i < m - 1; i++) {
+                for (int j = 1; j < n - 1; j++) {
+                    // calculate residual
+                    res.setDouble(i * n + j, (phi.getDouble((i + 1) * n + j) +
+                            phi.getDouble((i - 1) * n + j) - 2 * phi.getDouble(i * n + j)) /
+                            dy2.getDouble(i * n + j) + (phi.getDouble(i * n + j + 1) +
+                            phi.getDouble(i * n + j - 1) - 2 * phi.getDouble(i * n + j)) /
+                            dx2.getDouble(i * n + j) + array.getDouble(i * n + j));
+                    // iterator
+                    phi.setDouble(i * n + j, phi.getDouble(i * n + j) + (1 + sorIndex) * res.getDouble(i * n + j) /
+                            (2 / dx2.getDouble(i * n + j) + 2 / dy2.getDouble(i * n + j)));
+                }
+            }
+
+            if (ArrayMath.max(res).doubleValue() < epsilon) {
+                break;
+            }
+        }
+
+        return phi;
+    }
+
 }