add val to swm_cartesian.py

Author: cenamiller

swm_python/swm_cartesian.py

@@ -16,7 +16,8 @@
 L_OUT = True # args.L_OUT
 VIS = False
 VIS_DT=10
-ITMAX = 4000
+VAL=True
+ITMAX = 1000
 dt = 90.
 tdt = dt
 dx = 100000.
@@ -77,6 +78,32 @@ def live_plot3(fu, fv, fp, title=''):
     #plt.ylabel('y')
     plt.show()
 
+def read_arrays(u_file, v_file, p_file):
+    u_v = np.loadtxt(u_file,)
+    v_v = np.loadtxt(v_file)
+    p_v = np.loadtxt(p_file)
+    return u_v, v_v, p_v
+
+def live_plot_val(fu, fv, fp, title=''):
+    mxu = fu.max()
+    mxv = fv.max()
+    mxp = fp.max()
+    clear_output(wait=True)
+    fig, (ax1, ax2, ax3) = plt.subplots(figsize=(13, 3), ncols=3)
+
+    pos1 = ax1.imshow(fp, cmap='Blues', vmin=-mxp, vmax=mxp,interpolation='none')
+    ax1.set_title('p')
+    plt.colorbar(pos1,ax=ax1)
+    pos2 = ax2.imshow(fu, cmap='Reds', vmin=-mxu, vmax=mxu,interpolation='none')
+    ax2.set_title('u')
+    plt.colorbar(pos2,ax=ax2)
+    pos3 = ax3.imshow(fv, cmap='Greens',vmin=-mxv, vmax=mxv,interpolation='none')
+    ax3.set_title('v')
+    plt.colorbar(pos3, ax=ax3)
+
+    fig.suptitle(title)
+    plt.show()
+
 # Initial values of the stream function and p
 # for i in range(M + 1):
 #    for j in range(N + 1):
@@ -425,6 +452,34 @@ def copy_var(
     print("t200: ",dt2)
     print("t300: ",dt3)
 
+    if VAL:
+
+        u_val_f = 'ref/u.64.64.IT4000.txt'
+        v_val_f = 'ref/v.64.64.IT4000.txt'
+        p_val_f = 'ref/p.64.64.IT4000.txt'
+        uval = np.zeros((M_LEN, N_LEN))
+        vval = np.zeros((M_LEN, N_LEN))
+        pval = np.zeros((M_LEN, N_LEN))
+
+        uref, vref, pref = read_arrays(v_val_f, u_val_f, p_val_f)
+        uval = uref-unew[:,:,0]
+        vval = vref-vnew[:,:,0]
+        pval = pref-pnew[:,:,0]
+        
+        uLinfN= np.linalg.norm(uval, np.inf)
+        vLinfN= np.linalg.norm(vval, np.inf)
+        pLinfN= np.linalg.norm(pval, np.inf)
+
+        
+
+        live_plot_val(uval, vval, pval, "Val")
+        print("uLinfN: ", uLinfN)
+        print("vLinfN: ", vLinfN)
+        print("pLinfN: ", pLinfN)
+        print("udiff max: ",uval.max())
+        print("vdiff max: ",vval.max())
+        print("pdiff max: ",pval.max())
+
 
 # gt4py NEXT part!!!!!!!!!!!!!!!!!!!!
 
@@ -547,3 +602,31 @@ def calc_h_program(
            print(" diagonal elements of p:\n", pnew[:,:,0].diagonal()[:-1])
            print(" diagonal elements of u:\n", unew[:,:,0].diagonal()[:-1])
            print(" diagonal elements of v:\n", vnew[:,:,0].diagonal()[:-1])
+           
+    if VAL:
+        print("val")
+        u_val_f = 'ref/u.64.64.IT4000.txt'
+        v_val_f = 'ref/v.64.64.IT4000.txt'
+        p_val_f = 'ref/p.64.64.IT4000.txt'
+        uval = np.zeros((M_LEN, N_LEN))
+        vval = np.zeros((M_LEN, N_LEN))
+        pval = np.zeros((M_LEN, N_LEN))
+
+        uref, vref, pref = read_arrays(v_val_f, u_val_f, p_val_f)
+        uval = uref-unew
+        vval = vref-vnew
+        pval = pref-pnew
+        
+        uLinfN= np.linalg.norm(uval, np.inf)
+        vLinfN= np.linalg.norm(vval, np.inf)
+        pLinfN= np.linalg.norm(pval, np.inf)
+
+        
+
+        live_plot_val(uval, vval, pval, "Val")
+        print("uLinfN: ", uLinfN)
+        print("vLinfN: ", vLinfN)
+        print("pLinfN: ", pLinfN)
+        print("udiff max: ",uval.max())
+        print("vdiff max: ",vval.max())
+        print("pdiff max: ",pval.max())