convergence test for SW with no bathymetry no source no analytical solution

Test1D/Makefile.1D.gfortran

@@ -12,16 +12,31 @@ SRC = ../Src1D
 FFLAGS = -DLINUX -J$(MODDIR)  -cpp -c $(OPT) -ffree-line-length-none
 LDFLAGS= -J$(MODDIR)  -cpp $(OPT)  -ffree-line-length-none
 
-INIT= init_bc_euler
-MODEL_VAR= variable_def_euler
+#INIT= init_bc_euler
+#UTIL = utils
+#MODEL_VAR= variable_def_euler
 
 #INIT= init_bc_wave_1D
 #MODEL_VAR= variable_def_wave_1D
+#UTIL = utils
 
-#INIT= init_bc_scalar_1D
-#MODEL_VAR= variable_def_scalar_1D
+INIT= init_bc_sw
+MODEL_VAR= variable_def_sw
+UTIL = utils_sw
 
-OBJS = $(addprefix $(OBJDIR)/, elements_1D.o  param2d.o $(MODEL_VAR).o  aretes.o scheme.o overloading.o  Model.o geometry.o algebra.o utils.o postprocessing.o timestepping.o $(INIT).o precision.o)
+#INIT= init_bc_burgers
+#MODEL_VAR= variable_def_burgers
+#UTIL = utils
+
+#INIT= init_bc_scalar
+#MODEL_VAR= variable_def_scalar
+#UTIL = utils
+
+#INIT= init_bc_scalar
+#MODEL_VAR= variable_def_damped_scalar
+#UTIL = utils
+
+OBJS = $(addprefix $(OBJDIR)/, elements_1D.o  param2d.o $(MODEL_VAR).o  aretes.o scheme.o overloading.o  Model.o geometry.o algebra.o $(UTIL).o postprocessing.o timestepping.o $(INIT).o precision.o)
 
 
 dec: $(MODDIR) $(OBJDIR) $(BINDIR) $(OBJS) $(SRC)/main_dec.f90 
@@ -37,7 +52,7 @@ $(BINDIR):
 	mkdir -p $(BINDIR)	
 
 
-$(OBJDIR)/Model.o: $(SRC)/Model.f90 $(OBJDIR)/utils.o $(OBJDIR)/param2d.o $(OBJDIR)/precision.o
+$(OBJDIR)/Model.o: $(SRC)/Model.f90 $(OBJDIR)/param2d.o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/Model.o $(SRC)/Model.f90
 
 $(OBJDIR)/aretes.o: $(SRC)/aretes.f90 $(OBJDIR)/precision.o
@@ -49,7 +64,7 @@ $(OBJDIR)/overloading.o: $(SRC)/overloading.f90 $(OBJDIR)/$(MODEL_VAR).o $(OBJDI
 $(OBJDIR)/param2d.o: $(SRC)/param2d.f90 $(OBJDIR)/$(MODEL_VAR).o $(OBJDIR)/elements_1D.o $(OBJDIR)/aretes.o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/param2d.o $(SRC)/param2d.f90
 
-$(OBJDIR)/geometry.o: $(SRC)/param2d.f90 $(OBJDIR)/elements_1D.o $(SRC)/geometry.f90 $(OBJDIR)/precision.o
+$(OBJDIR)/geometry.o: $(SRC)/param2d.f90 $(OBJDIR)/elements_1D.o $(OBJDIR)/$(INIT).o $(SRC)/geometry.f90 $(OBJDIR)/precision.o
 	 $(F90) $(FFLAGS) -o $(OBJDIR)/geometry.o $(SRC)/geometry.f90
 
 $(OBJDIR)/algebra.o: $(SRC)/algebra.f90 $(OBJDIR)/precision.o
@@ -58,22 +73,22 @@ $(OBJDIR)/algebra.o: $(SRC)/algebra.f90 $(OBJDIR)/precision.o
 $(OBJDIR)/elements_1D.o: $(SRC)/elements_1D.f90 $(OBJDIR)/$(MODEL_VAR).o $(OBJDIR)/algebra.o $(OBJDIR)/overloading.o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/elements_1D.o $(SRC)/elements_1D.f90
 
-$(OBJDIR)/$(MODEL_VAR).o: $(OBJDIR)/algebra.o $(SRC)/$(MODEL_VAR).f90   $(OBJDIR)/precision.o
+$(OBJDIR)/$(MODEL_VAR).o: $(OBJDIR)/algebra.o $(SRC)/$(MODEL_VAR).f90   $(OBJDIR)/$(UTIL).o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/$(MODEL_VAR).o $(SRC)/$(MODEL_VAR).f90
 
 $(OBJDIR)/scheme.o: $(SRC)/scheme.f90  $(OBJDIR)/$(MODEL_VAR).o $(OBJDIR)/elements_1D.o $(OBJDIR)/overloading.o $(OBJDIR)/aretes.o $(OBJDIR)/Model.o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/scheme.o $(SRC)/scheme.f90
 
-$(OBJDIR)/utils.o:$(SRC)/utils.f90 $(OBJDIR)/elements_1D.o $(OBJDIR)/$(MODEL_VAR).o $(OBJDIR)/$(INIT).o  $(OBJDIR)/precision.o
-	$(F90) $(FFLAGS) -o $(OBJDIR)/utils.o $(SRC)/utils.f90
+$(OBJDIR)/$(UTIL).o:$(SRC)/$(UTIL).f90 $(OBJDIR)/precision.o
+	$(F90) $(FFLAGS) -o $(OBJDIR)/$(UTIL).o $(SRC)/$(UTIL).f90
 
-$(OBJDIR)/$(INIT).o: $(SRC)/$(INIT).f90 $(OBJDIR)/param2d.o $(OBJDIR)/overloading.o $(OBJDIR)/precision.o
+$(OBJDIR)/$(INIT).o: $(SRC)/$(INIT).f90 $(OBJDIR)/param2d.o $(OBJDIR)/overloading.o $(OBJDIR)/$(UTIL).o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/$(INIT).o $(SRC)/$(INIT).f90
 
 $(OBJDIR)/timestepping.o: $(SRC)/timestepping.f90 $(OBJDIR)/overloading.o $(OBJDIR)/elements_1D.o $(OBJDIR)/$(MODEL_VAR).o $(OBJDIR)/param2d.o $(OBJDIR)/scheme.o $(OBJDIR)/Model.o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/timestepping.o $(SRC)/timestepping.f90
 
-$(OBJDIR)/postprocessing.o: $(SRC)/postprocessing.f90 $(OBJDIR)/param2d.o $(OBJDIR)/utils.o $(OBJDIR)/Model.o $(OBJDIR)/precision.o
+$(OBJDIR)/postprocessing.o: $(SRC)/postprocessing.f90 $(OBJDIR)/param2d.o $(OBJDIR)/$(INIT).o $(OBJDIR)/Model.o $(OBJDIR)/precision.o
 	$(F90) $(FFLAGS) -o $(OBJDIR)/postprocessing.o $(SRC)/postprocessing.f90
 	
 $(OBJDIR)/precision.o: $(SRC)/precision.f90
@@ -86,6 +101,3 @@ clean:
 	rm $(SRC)/*.f90~ 
 	rm *.mod
 
-
-
-

Test1D/Test1D_SW_convergence_no_bath/Data/don1d

+1000   nt
+12    itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
+3 3  ordre 3: # of iteration for DEC 
+4   scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
+0.119   theta parameter in Burman stabilization term 
+0.2   theta2 parameter in Burman stabilization second derivative term 
+0.1    cfl 
+100000    ktmax
+3.0     tmax
+100     ifre
+42       test case SW
+.FALSE.   !MOOD

Test1D/Test1D_SW_convergence_no_bath/convergence_no_an_Pgl.py

+#!/usr/bin/python
+import numpy as np 
+from numpy import linalg as LA
+import scipy
+import matplotlib.pyplot as pl
+import sys 
+import itertools
+#import exact_sol as ex
+
+# run the convergence without an exact solution
+# files saved in columns x, h, u
+# ns and the file are below and must be changed accordingly
+
+def remove_extra(x, u):
+   ind=[]
+   for i in range(len(x)-1):
+      if x[i]==x[i+1]:
+          ind.append(i)
+   x=np.delete(x,ind,0)
+   u=np.delete(u,ind,0)
+   return x,u
+
+pl.rc("text", usetex=True)
+pl.rc("font", family="serif")
+
+# Grids for convergence study
+#nNodes = [5,10,20,40,80,160,320,640,1280,2560]
+#nNodes = [5,10,20,40,80,160,320,640,1280]
+#nNodes = [5,10,20,40,80,160,320,640]
+#nNodes = [5,10,20,40,80,160,320]
+#nNodes = [5,10,20,40,80,160]
+#nNodes = [5,10,20,40,80]
+#nNodes = [5,10,20,40]
+nNodes = [8,16,32,64,128]#,256,512, 1024]#,2048,4096]
+#nNodes = [25,50,100,200]#,400,800]
+# Test name
+test = "test"
+
+# scheme
+scheme = "scheme4"
+
+# Scheme order
+orders = [1, 2,3]#,4]
+
+marker = itertools.cycle(('s','o','*'))
+
+varnames = ['h','u']
+
+for nvar, var in enumerate(varnames):
+		 # Initialize plots
+		 fig_err_L1 = pl.figure()
+		 ax_err_L1  = fig_err_L1.add_subplot(1,1,1)
+		
+		 # Initialize array to store errors
+		 err_L1 = np.zeros((np.size(orders),np.size(nNodes)-1))
+		 dx     = np.zeros((np.size(orders),np.size(nNodes)-1))
+
+		 # Initialize array to store orders
+		 order_L1 = np.zeros((np.size(orders),np.size(nNodes)-1))
+		
+		 for k, order in enumerate(orders):
+		    xs =[]
+		    num_sols = []
+		    for i, N in enumerate(nNodes):
+
+		       # Load numerical solution
+		       fname = "Pgl"+str(order)+"/"+str(test)+"_"+scheme+"_"+str(N)+".dat"  #here the folder and filenames
+		       x, num_sol = np.loadtxt(fname,delimiter=None,usecols=(0,nvar+1), unpack=True)
+		       x, num_sol = remove_extra(x,num_sol)
+		       xs.append(x)
+		       num_sols.append(num_sol)
+
+		    for i in range(len(nNodes)-1):
+		       err_L1[k,i]=LA.norm(num_sols[i][::order]-num_sols[i+1][::2*order],ord=2)/float(np.sqrt(nNodes[i]))
+
+		       # Compute exact solution
+	#         ex_sol = ex.exact_sw_smooth(x)
+	 
+		       # Compute errors
+		       #err_L1[k,i] = LA.norm(num_sol-ex_sol[nvar],ord=np.inf)
+	#         err_L1[k,i] = LA.norm(num_sol-ex_sol[nvar],ord=1)/float(N)
+		       #err_L1[k,i] = LA.norm(num_sol-ex_sol[nvar],ord=2)/float(np.sqrt(N))
+		       #err_L1[k,i] = LA.norm((num_sol-ex_sol[nvar])*dxvec,ord=1)
+		       
+		       if (i > 0):
+		          order_L1[k,i] = -( np.log(err_L1[k][i])-np.log(err_L1[k][i-1]) ) / ( np.log(nNodes[i])-np.log(nNodes[i-1]) )
+		    
+		    print("Variable : ", var)
+		    print("L1-errors for order ", order, " :")
+		    print("Computed errors : ")
+		    print(err_L1[k,:])
+		    print("Computed orders : ")
+		    print(order_L1[k,:])
+
+		    # Plot figures
+
+		    # error vs mesh size
+		    ax_err_L1.loglog(nNodes[1:],err_L1[k,:],"-"+marker.next(),label="Pgl"+str(order))
+		    ref_err = [err_L1[k,1]*nNodes[1]**(float(order+1))*nNodes[i]**(-float(order+1)) for i in range(np.size(nNodes[1:]))]
+		    ax_err_L1.loglog(nNodes[1:],ref_err,"--",label="order "+str(order+1))
+		    ax_err_L1.set_xlabel(r"$N$")
+		    ax_err_L1.set_ylabel(r"$L_2$-error")
+		    pl.xlim(nNodes[1]*0.7,nNodes[-1]*1.3)
+		    ax_err_L1.legend(loc=0)
+		    pl.title("Convergence of "+var)      
+
+     # Save plots
+		 fig_err_L1.savefig(scheme+"_"+var+"_convergence_no_an_Pgl.pdf",format="pdf")
+
+pl.show()

Test1D/Test1D_SW_convergence_no_bath/run_tests_Pgl.py

+import numpy as np
+import shutil, os
+import re
+from joblib import Parallel, delayed
+import multiprocessing
+
+#In this function we run the tests for all ns for one parameter (e.g. B2)
+
+def run_single_test(cand, param, test_folder):
+	instruction ="" 
+	foldName = "Pgl"+str(param)
+	instruction +="mkdir "+foldName+" \n"
+	instruction +="mkdir "+foldName+"/Data \n"
+	instruction +="cp Data/don1d "+foldName+"/Data/don1d \n"
+	os.system(instruction)
+	modify_don(foldName+"/Data/don1d", param)
+	ns= [int(k) for k in 2.**np.linspace(3.,10.,8)]  #here are the ns
+	for n in ns:
+#		n=param[1]
+		instruction ="" 
+		instruction +="cd "+foldName+" \n "  #changing folder into Bn
+		instruction += "pwd \n"
+		instruction +=  "../../bin1D/main_dec.out "+str(n)+ " \n"  #running the test with n. call as if we are in the Bn folder
+		instruction += "mv TEST/sol_prim_last.dat test_scheme4_"+str(n)+".dat  \n"
+		instruction += "mv TEST/sol_cons_last.dat test_cons_scheme4_"+str(n)+".dat  \n"
+		instruction += "cd TEST  \n"
+		instruction += "rm -rf * \n"
+		os.system(instruction)
+	return
+
+
+#This is the function to run more paramaters
+def modify_don(fname, param):
+	don=[]
+	with open(fname, 'r') as a:
+		for idx, line in enumerate(a.readlines()):
+			don.append(line)
+
+	types=[11,12,13,14]
+	order = [[2,2], [3,3], [4,4], [5,5]]
+	burman=[[1., 0.],[0.1,0.1], [0.1,0.1],[0.1,0.1]]
+	burman=[[0.119, 0],[0.00346,0],[0.000113,0],[0.1,0 ]]#[[0.05623, 0.],[0.0056,  0], [0.001778, 0],[0.1,0 ]]
+	burman=[[0.119, 0],[0.05,0],[0.05,0],[0.01,0 ]]#[[0.05623, 0.],[0.0056,  0], [0.001778, 0],[0.1,0 ]]
+	#burman=[[0.119, 0],[0.0,0],[0.0,0],[0.0,0 ]]#[[0.05623, 0.],[0.0056,  0], [0.001778, 0],[0.1,0 ]]
+	CFLs=[0.6, 0.2, 0.1, 0.05]
+
+	pol=int(param)-1
+	f=open(fname,'w')
+	for idx, line in enumerate(don):
+		if idx==2:
+			f.write(str(order[pol][0]) +' '+str(order[pol][1])+'  ordre 3: # of iteration for DEC \n')
+		elif idx==1:				
+			f.write(str(types[pol])+'    itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3\n')
+		elif idx==4:				
+			f.write(str(burman[pol][0])+'   theta parameter in Burman stabilization term \n')
+		elif idx==5:				
+			f.write(str(burman[pol][1])+'   theta2 parameter in Burman stabilization second derivative term \n')
+		elif idx==6:				
+			f.write(str(CFLs[pol])+'   cfl \n')
+		else:
+			f.write(don[idx])
+	f.close()
+
+	return
+
+
+
+
+#vary B1/B2/B3
+#vary theta burman between 0, 0.5
+#vary alpha lax friedrics
+test_folder = ""
+
+num_cores = min(multiprocessing.cpu_count(), 4)
+
+param_domain=[1,2,3,4] #[1,1]]#,[2,1],[3,1],[4,1]]#, [1]]#,[2]] [1,0],[2,0],[3,0],[4,0],
+
+Parallel(n_jobs=num_cores)(delayed(run_single_test)(cand, param_domain[cand], test_folder) for cand in range(len(param_domain)))
+
+
+
+