added error analysis for vortex case

mhd/error_analysis.py

+import numpy as np
+import argparse
+from lbm_cl import *
+from default_variables import *
+from matplotlib import pyplot as plt
+
+
+m = 9
+
+
+
+def L1_error (error_field, macrov, nx, ny):
+
+    reshaped_data = np.reshape(error_field, (4, m, 
+                               nx, ny))
+    macrov_data = np.sum(error_field[:, macrov, :, :], axis=0)
+    plt.imshow(macrov_data, interpolation='bilinear')
+    plt.colorbar()
+    plt.show()
+
+    l1err = 0.0
+
+    for wx in macrov_data:
+    	for w in wx:
+    		#print(np.abs(w))
+    		l1err += np.abs(w);
+
+    return l1err / (nx*ny)
+
+
+def L2_error (error_field, macrov, nx, ny):
+
+    reshaped_data = np.reshape(error_field, (4, m, 
+                               nx, ny))
+    macrov_data = np.sum(error_field[:, macrov, :, :], axis=0)
+
+    l2err = 0.0
+
+    for wx in macrov_data:
+    	for w in wx:
+    		#print(np.abs(w))
+    		l2err += w**2;
+
+    return np.sqrt(l2err) / np.sqrt(100*nx*ny)
+
+
+
+
+if __name__ == '__main__':
+	print ("analysis of the vortex case !")
+
+
+	parser = argparse.ArgumentParser(description="perform error analysis")
+	parser.add_argument('-nm', '--nmini', type=int, default=128, 
+		                help='minimum resolution for convergence')
+	parser.add_argument('-nt', '--ntot', type=int, default=1, help='number of grids for convergence')
+	args = parser.parse_args()
+	nx    = args.nmini
+	nconv = args.ntot 
+
+	#nconv = 4
+	errors=np.zeros(nconv)
+	#nx = 128
+
+	print(nconv, nx)
+
+	for i in range (4):
+		ny = nx
+		x, y, wplot_gpu, error, tcpu = solve_ocl (nx=nx, ny=ny, Lx=10, Ly=10, 
+			Tmax=1.0, test="vortex", savemacrodata=False, savekineticdata=False,
+			computeerror=True)
+
+		errors[i] = L2_error (error, eBx, nx, ny)
+
+		print ("error results:")
+		print ("(",nx,",",ny,") : ",errors[i])
+
+		nx *=2
+
+	for errs in errors:
+		print(errs)

mhd/lbm_cl.py

@@ -64,7 +64,7 @@ def set_parameters (m, n, nx, ny, Lx, Ly):
 
 def solve_ocl(m=m_, n=n_, nx=nx_, ny=ny_, Lx=Lx_, Ly=Ly_, Tmax=Tmax_, test=test_,
               animate=False, precision="single",savekineticdata="False", 
-              savemacrodata="False"):
+              savemacrodata="False", computeerror=False):
 
 
     testname = test
@@ -98,6 +98,8 @@ def solve_ocl(m=m_, n=n_, nx=nx_, ny=ny_, Lx=Lx_, Ly=Ly_, Tmax=Tmax_, test=test_
     fn_gpu      = cl.Buffer(ctx, mf.READ_WRITE, size=buffer_size)
     fnp1_gpu    = cl.Buffer(ctx, mf.READ_WRITE, size=buffer_size)
     kinetic_gpu = cl.Buffer(ctx, mf.READ_WRITE, size=kinetic_buffer_size)
+    if computeerror:
+        error_gpu   = cl.Buffer(ctx, mf.READ_WRITE, size=buffer_size) 
     #divb_gpu = cl.Buffer(ctx, mf.READ_WRITE, size=buffer_size)
                          
     if (nx>parameters.get('ncpr')):
@@ -121,7 +123,12 @@ def solve_ocl(m=m_, n=n_, nx=nx_, ny=ny_, Lx=Lx_, Ly=Ly_, Tmax=Tmax_, test=test_
     elapsed = 0.
     
     fn_cpu      = np.zeros((4 * m * ncprx * ncpry, ), dtype=np_real)
+    if computeerror:
+        error_cpu   = np.zeros((4 * m * nx * ny, ), dtype=np_real) 
+    else:
+        error_cpu = 0.0
     kinetic_cpu = np.zeros((1 * ncprx * ncpry, ), dtype=np_real)
+    time_vec    = np.zeros((1,), dtype=np_real)
     #divb_cpu = np.zeros((1 * nx * ny, ), dtype=np_real)
 
     if animate:
@@ -163,6 +170,7 @@ def solve_ocl(m=m_, n=n_, nx=nx_, ny=ny_, Lx=Lx_, Ly=Ly_, Tmax=Tmax_, test=test_
                     pstmgr.h5_savedata(fn_cpu, kinetic_cpu, parameters, t, ite)
                 
         t += parameters.get('dt')
+        time_vec[0] = t
         #event = prg.time_step(queue, (nx * ny, ), None, wn_gpu, wnp1_gpu)
         event = prg.time_step(queue, (nx * ny, ), None, fn_gpu, fnp1_gpu)
         #event = prg.time_step(queue, (nx * ny, ), None, wn_gpu, wnp1_gpu, wait_for = [event])
@@ -180,21 +188,30 @@ def solve_ocl(m=m_, n=n_, nx=nx_, ny=ny_, Lx=Lx_, Ly=Ly_, Tmax=Tmax_, test=test_
         ite += 1
 
     # copy OpenCL data to CPU and return the results
+
+    if computeerror:
+        event = prg.error_field (queue, (nx * ny,), None, fnp1_gpu, error_gpu, time_vec[0]);
+        cl.enqueue_copy(queue, error_cpu, error_gpu).wait()
+        queue.finish()
+        error_cpu = np.reshape (error_cpu, (4,m,nx,ny))
+
     cl.enqueue_copy(queue, fn_cpu, fnp1_gpu).wait()
     queue.finish()
 
     wplot_gpu = np.reshape(fn_cpu, (4, m, ncprx, ncpry))
     plt.show(block=True)
-    
-    return pstmgr.x, pstmgr.y, wplot_gpu
+
+
+    return pstmgr.x, pstmgr.y, wplot_gpu, error_cpu, elapsed
 
 
 if __name__ == '__main__':
 
     args = parse_cl_args(n=nx_, L=Lx_, tmax=Tmax_, test=test_,
                          help_test = help_test_,
-                         description='Solve orszag-tang using LBM on PyOpenCL')
+                         description='Solve orszag-tang using LBM on PyOpenCL',
+                         computeerror = False)
 
     # gpu solve
-    x, y, wplot_gpu = solve_ocl(**vars(args))
+    x, y, wplot_gpu, error, tcpu = solve_ocl(**vars(args))
 

mhd/lbm_common.h

@@ -160,7 +160,7 @@ void conservatives(real *y, real *w) {
 
   w[eRho] = y[eRho];   // rho
   w[eU] = y[eRho]*y[eU]; // rho u
-  w[eE] = y[eE]/(gam-1) + y[eRho]*(y[eU]*y[eU]+y[eV]*y[eV]+y[eW]*y[eW])/2
+  w[eE] = y[eP]/(gam-1) + y[eRho]*(y[eU]*y[eU]+y[eV]*y[eV]+y[eW]*y[eW])/2
     + (y[eBx]*y[eBx]+y[eBy]*y[eBy]+y[eBz]*y[eBz])/2; // rho E
   w[eV] = y[eRho]*y[eV];  // rho v
   w[eW] = y[eRho]*y[eW]; // rho w

mhd/vortex_kernels.cl

@@ -18,7 +18,7 @@
 
 #include "lbm_common.h"
 
-void exact_sol(real* x, real t, real* w){
+void exact_sol(real* x, const real t, real* w){
 
   real vx = 1.0_F; // vortex displacement velocity
   real vy = 1.0_F;
@@ -47,7 +47,7 @@ void exact_sol(real* x, real t, real* w){
   real du = dvelo*(yt-x[eY]);
   real dv = dvelo*(x[eX]-xt);
   
-  real dB = mu/(2.0*M_PI);
+  real dB = mu/(4.0*M_PI*sqrt(M_PI));
   dB *= exp(q*(1.0-r2));
   real dBx = dB*(yt-x[eY]);
   real dBy = dB*(x[eX]-xt);
@@ -85,7 +85,7 @@ __kernel void init_sol(__global  real *fn){
 
   real wnow[_M];
 
-  real t = 2.5;
+  const real t = 0.0;
   real xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
   
   exact_sol(xy, t, wnow);
@@ -100,4 +100,32 @@ __kernel void init_sol(__global  real *fn){
     fn[imem] =  fnow[ik];
     //fn[imem] = j;
   }
+}
+
+
+// error computation
+__kernel void error_field (__global const real *fn, __global real *dfn, const real t)
+{
+  int id = get_global_id(0);
+  
+  int i = id % _NX;
+  int j = id / _NX;
+
+  int ngrid = _NX * _NY;
+
+  real wref[_M];
+  real xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
+  
+  exact_sol(xy, t, wref);
+  
+  real fref[_N];
+  w2f(wref, fref);
+
+  //printf("x=%f, y=%f \n",xy[0],xy[1]);
+  // load middle value
+  for(int ik = 0; ik < _N; ik++){
+    int imem = i + j * _NX + ik * ngrid;
+    dfn[imem] =  fn[imem]-fref[ik];
+    //fn[imem] = j;
+  }
 }
\ No newline at end of file

utils/utils.py

@@ -180,6 +180,8 @@ def parse_cl_args(n=256, L=1.0, tmax=1.0, test='None',
                         help = 'floating point precision')
     parser.add_argument('-test', '--test', type=str, default=test,
                         help = help_test)
+    parser.add_argument('-er', '--error', action='store_true',
+                        help='compute the solution error field')
     args = parser.parse_args()
 
     # Prepare to build a square domain with nx = ny = n