Commit ba819485 authored by Paola Bacigaluppi's avatar Paola Bacigaluppi
Browse files

Updated the Readme, The test 1D Sod problem with some files (run_tests) to run...

Updated the Readme, The test 1D Sod problem with some files (run_tests) to run the test on diverse mesh sizes by doing e.g ./run_tests Sod B1 from the B1 folder, and a python script to plot B1,B2,B3 on the diverse meshes and for different variables
parent f9d52e4d
100 nt
100 nt
1 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
2 2 ordre 3: # of iteration for DEC
5 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
1. alpha_jump: parameter used in Burman stabilization term
0.0 alpha_jump: parameter used in Burman stabilization term
4 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
1. alpha_jump: parameter used in Burman stabilization term
0. alpha_jump: parameter used in Burman stabilization term
0.1 cfl
100000 ktmax
3.0 tmax
......
#! /bin/bash
test=$1
basis=$2
for mesh in 50 100 200 400
do
./../../bin1D/main_dec.out $mesh
mv ./TEST/sol_prim_last.dat ./TEST/test_$test"_"$basis"_Mood_"$mesh".dat"
mv ./TEST/FluxIndicator_last.dat ./TEST/test_$test"_"$basis"_FluxIndicator_"$mesh".dat"
mv ./TEST/DiagIndicator_last.dat ./TEST/test_$test"_"$basis"_DiagIndicator_"$mesh".dat"
rm -rf ./TEST/sol_prim*
rm -rf ./TEST/sol_cons*
rm -rf ./TEST/Flux*
rm -rf ./TEST/Entrop*
rm -rf ./TEST/Diag*
done
#! /bin/bash
test=$1
basis=$2
for mesh in 50 100 200 400
do
./../../bin1D/main_dec.out $mesh
mv ./TEST/sol_prim_last.dat ./TEST/test$test"_"$basis"_Mood_"$mesh".dat"
mv ./TEST/FluxIndicator_last.dat ./TEST/test$test"_"$basis"_FluxIndicator_"$mesh".dat"
mv ./TEST/DiagIndicator_last.dat ./TEST/test$test"_"$basis"_DiagIndicator_"$mesh".dat"
rm -rf ./TEST/sol_prim*
rm -rf ./TEST/sol_cons*
rm -rf ./TEST/Flux*
rm -rf ./TEST/Entrop*
rm -rf ./TEST/Diag*
done
import numpy as np
import shutil, os
import re
from joblib import Parallel, delayed
import multiprocessing
def run_single_test(cand, param):
instruction = "mkdir param"+str(cand)+" \n "
instruction +="cd param"+str(cand)+" \n "
instruction +="mkdir Data"+" \n "
instruction +="cp ../Data/don1d Data/don1d"+" \n "
#instruction +="cp ../../Data/parameter Data/parameter"+" \n "
#mesh_name = "R13_o"+str(int(param[0]))+param[2]
#instruction +="cp ../../"+mesh_name+".msh mesh.msh"+" \n "
os.system(instruction)
modify_don( "param"+str(cand)+"/Data/don1d", param)
instruction=""
instruction+="cd param"+str(cand)+"/" +" \n "
instruction += "../../../bin1D_pb/main_dec.out"
os.system(instruction)
return
def modify_don(fname, param):
don=[]
with open(fname, 'r') as a:
for idx, line in enumerate(a.readlines()):
don.append(line)
types=[1,2,5]
order = [[2,3], [3,4], [4,5]]
pol=int(param[0])
f=open(fname,'w')
for idx, line in enumerate(don):
if idx==1:
f.write(str(types[pol-1]) + " Itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3 \n")
elif idx==2:
f.write(str(order[pol-1][0])+' '+str(order[pol-1][1])+' ordre 3: # of iteration for DEC \n')
elif idx==4:
f.write(str(param[1])+' alphajump parameter in Burman stabilization term \n')
elif idx==5:
f.write(str(param[2])+' alphajump2 parameter in Burman stabilization term \n')
else:
f.write(don[idx])
f.close()
return
#vary B1/B2/B3
#vary theta burman between 0, 0.5
#vary alpha lax friedrics
order=2 #B1 is 1 B2 is 2 B3is 3
num_cores = min(multiprocessing.cpu_count(), 60)
interval_min = [ 0.0 , 0.0]
interval_max = [ 1., 1.]
num_par = [ 6 , 6]
grids=[]
for j in range(2):
grids.append(np.linspace(interval_min[j],interval_max[j],num_par[j]))
param_domain=[]
for b in range(num_par[0]):
for d in range(num_par[1]):
param_domain.append([order,grids[0][b],grids[1][d]])
Parallel(n_jobs=num_cores)(delayed(run_single_test)(cand, param_domain[cand]) for cand in range(len(param_domain)))
100 nt
3 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
2 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
3 3 ordre 3: # of iteration for DEC
5 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
4 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
1. alpha_jump: parameter used in Burman stabilization term
0.0 alpha_jump: parameter used in Burman stabilization term
0.1 cfl
100000 ktmax
100000 ktmax
3.0 tmax
100 ifre
1 test case Euler: 0 - smooth (convergence), 1 - Sod, 2 - Shu-Osher, 3 - Woodward-Colella, 6 - Woodward-Colella left
.true.
2
5
-1
100 nt
2 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
3 3 ordre 3: # of iteration for DEC
5 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
4 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
1. alpha_jump: parameter used in Burman stabilization term
0.0 alpha_jump: parameter used in Burman stabilization term
0.1 cfl
100000 ktmax
100000 ktmax
3.0 tmax
100 ifre
1 test case Euler: 0 - smooth (convergence), 1 - Sod, 2 - Shu-Osher, 3 - Woodward-Colella, 6 - Woodward-Colella left
......
import numpy as np
import re
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
def colors(theta):
return (theta[1], 0.5,theta[2])
order=2
interval_min = [ 0.0 , 0.0]
interval_max = [ 1., 1.]
num_par = [ 6 , 6]
grids=[]
for j in range(2):
grids.append(np.linspace(interval_min[j],interval_max[j],num_par[j]))
param_domain=[]
for b in range(num_par[0]):
for d in range(num_par[1]):
param_domain.append([order,grids[0][b],grids[1][d]])
#selection=range(5,len(param_domain),6)
#selection=range(30,36)
selection=range(6,36)
for cand in selection:
par=param_domain[cand]
solut=np.loadtxt("param"+str(cand)+"/sol_prim_last.dat")
par=param_domain[cand]
plt.plot(solut[:,0],solut[:,1], label="alpha1="+str(par[1])+"alpha2="+str(par[2]),color=colors(par))
solut=np.loadtxt("test1_exact_rho.dat")
plt.plot(solut[:,0],solut[:,1],label="exact", color=[0,0,0])
plt.title("Different alphas")
plt.legend()
plt.show()
plt.savefig("alpha_comparison_B"+str(par[0])+".pdf")
import numpy as np
import shutil, os
import re
from joblib import Parallel, delayed
import multiprocessing
def run_single_test(cand, param):
instruction = "mkdir param"+str(cand)+" \n "
instruction +="cd param"+str(cand)+" \n "
instruction +="mkdir Data"+" \n "
instruction +="cp ../Data/don1d Data/don1d"+" \n "
#instruction +="cp ../../Data/parameter Data/parameter"+" \n "
#mesh_name = "R13_o"+str(int(param[0]))+param[2]
#instruction +="cp ../../"+mesh_name+".msh mesh.msh"+" \n "
os.system(instruction)
modify_don( "param"+str(cand)+"/Data/don1d", param)
instruction=""
instruction+="cd param"+str(cand)+"/" +" \n "
instruction += "../../../bin1D_pb/main_dec.out"
os.system(instruction)
return
def modify_don(fname, param):
don=[]
with open(fname, 'r') as a:
for idx, line in enumerate(a.readlines()):
don.append(line)
types=[1,2,5]
order = [[2,3], [3,4], [4,5]]
pol=int(param[0])
f=open(fname,'w')
for idx, line in enumerate(don):
if idx==1:
f.write(str(types[pol-1]) + " Itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3 \n")
elif idx==2:
f.write(str(order[pol-1][0])+' '+str(order[pol-1][1])+' ordre 3: # of iteration for DEC \n')
elif idx==4:
f.write(str(param[1])+' alphajump parameter in Burman stabilization term \n')
elif idx==5:
f.write(str(param[2])+' alphajump2 parameter in Burman stabilization term \n')
else:
f.write(don[idx])
f.close()
return
#vary B1/B2/B3
#vary theta burman between 0, 0.5
#vary alpha lax friedrics
order=2 #B1 is 1 B2 is 2 B3is 3
num_cores = min(multiprocessing.cpu_count(), 60)
interval_min = [ 0.0 , 0.0]
interval_max = [ 1., 1.]
num_par = [ 6 , 6]
grids=[]
for j in range(2):
grids.append(np.linspace(interval_min[j],interval_max[j],num_par[j]))
param_domain=[]
for b in range(num_par[0]):
for d in range(num_par[1]):
param_domain.append([order,grids[0][b],grids[1][d]])
Parallel(n_jobs=num_cores)(delayed(run_single_test)(cand, param_domain[cand]) for cand in range(len(param_domain)))
#! /bin/bash
test=$1
basis=$2
for mesh in 50 100 200 400
do
./../../bin1D/main_dec.out $mesh
mv ./TEST/sol_prim_last.dat ./TEST/test_$test"_"$basis"_Mood_"$mesh".dat"
mv ./TEST/FluxIndicator_last.dat ./TEST/test_$test"_"$basis"_FluxIndicator_"$mesh".dat"
mv ./TEST/DiagIndicator_last.dat ./TEST/test_$test"_"$basis"_DiagIndicator_"$mesh".dat"
rm -rf ./TEST/sol_prim*
rm -rf ./TEST/sol_cons*
rm -rf ./TEST/Flux*
rm -rf ./TEST/Entrop*
rm -rf ./TEST/Diag*
done
#! /bin/bash
test=$1
basis=$2
for mesh in 50 100 200 400
do
./../../bin1D/main_dec.out $mesh
mv ./TEST/sol_prim_last.dat ./TEST/test$test"_"$basis"_Mood_"$mesh".dat"
mv ./TEST/FluxIndicator_last.dat ./TEST/test$test"_"$basis"_FluxIndicator_"$mesh".dat"
mv ./TEST/DiagIndicator_last.dat ./TEST/test$test"_"$basis"_DiagIndicator_"$mesh".dat"
rm -rf ./TEST/sol_prim*
rm -rf ./TEST/sol_cons*
rm -rf ./TEST/Flux*
rm -rf ./TEST/Entrop*
rm -rf ./TEST/Diag*
done
100 nt
5 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
5 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
4 4 ordre 3: # of iteration for DEC
5 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
1. alpha_jump: parameter used in Burman stabilization term
4 scheme: 1=supg, 2=psi, 3=mix, 4: galerkin+jump, 5: psi+jump 6 blend+jump 7: psi+galerkin2??
3. alpha_jump: parameter used in Burman stabilization term
10. alpha_jump: parameter used in Burman stabilization term
0.1 cfl
20000 ktmax
10000 ktmax
3.0 tmax
100 ifre
1 test case Euler: 0 - smooth (convergence), 1 - Sod, 2 - Shu-Osher, 3 - Woodward-Colella, 6 - Woodward-Colella left
.true.
2
5
-1
100 nt
5 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
4 4 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??
3. alpha_jump: parameter used in Burman stabilization term
10. alpha_jump: parameter used in Burman stabilization term
0.1 cfl
11 ktmax
3.0 tmax
1 ifre
1 test case Euler: 0 - smooth (convergence), 1 - Sod, 2 - Shu-Osher, 3 - Woodward-Colella, 6 - Woodward-Colella left
.true.
2
5
-1
import numpy as np
import re
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
#def colors(theta):
# return (theta[1], 0.5*theta[2],theta[2])
order=2
interval_min = [ 0.0 , 0.0]
interval_max = [ 1., 1.]
num_par = [ 6 , 6]
grids=[]
for j in range(2):
grids.append(np.linspace(interval_min[j],interval_max[j],num_par[j]))
param_domain=[]
for b in range(num_par[0]):
for d in range(num_par[1]):
param_domain.append([order,grids[0][b],grids[1][d]])
#selection=range(5,len(param_domain),6)
#selection=range(30,36)
selection=range(0,36)
for cand in selection:
par=param_domain[cand]
solut=np.loadtxt("param"+str(cand)+"/sol_prim_last.dat")
par=param_domain[cand]
plt.plot(solut[:,0],solut[:,1], label="alpha1="+str(par[1])+"alpha2="+str(par[2]))#,co lor=colors(par))
solut=np.loadtxt("test1_exact_rho.dat")
plt.plot(solut[:,0],solut[:,1],label="exact", color=[0,0,0])
plt.title("Different alphas")
plt.legend()
plt.show()
plt.savefig("alpha_comparison_B"+str(par[0])+".pdf")
import numpy as np
import shutil, os
import re
from joblib import Parallel, delayed
import multiprocessing
def run_single_test(cand, param):
instruction = "mkdir param"+str(cand)+" \n "
instruction +="cd param"+str(cand)+" \n "
instruction +="mkdir Data"+" \n "
instruction +="cp ../Data/don1d Data/don1d"+" \n "
#instruction +="cp ../../Data/parameter Data/parameter"+" \n "
#mesh_name = "R13_o"+str(int(param[0]))+param[2]
#instruction +="cp ../../"+mesh_name+".msh mesh.msh"+" \n "
os.system(instruction)
modify_don( "param"+str(cand)+"/Data/don1d", param)
instruction=""
instruction+="cd param"+str(cand)+"/" +" \n "
instruction += "../../../bin1D/main_dec.out"
os.system(instruction)
return
def modify_don(fname, param):
don=[]
with open(fname, 'r') as a:
for idx, line in enumerate(a.readlines()):
don.append(line)
types=[1,2,5]
order = [[2,3], [3,4], [4,5]]
pol=int(param[0])
f=open(fname,'w')
for idx, line in enumerate(don):
if idx==1:
f.write(str(types[pol-1]) + " Itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3 \n")
elif idx==2:
f.write(str(order[pol-1][0])+' '+str(order[pol-1][1])+' ordre 3: # of iteration for DEC \n')
elif idx==4:
f.write(str(param[1])+' alphajump parameter in Burman stabilization term \n')
elif idx==5:
f.write(str(param[2])+' alphajump2 parameter in Burman stabilization term \n')
else:
f.write(don[idx])
f.close()
return
#vary B1/B2/B3
#vary theta burman between 0, 0.5
#vary alpha lax friedrics
order=3 #B1 is 1 B2 is 2 B3is 3
num_cores = min(multiprocessing.cpu_count(), 60)
interval_min = [ 0.0 , 0.0]
interval_max = [ 1., 1.]
num_par = [ 6 , 6]
grids=[]
for j in range(2):
grids.append(np.linspace(interval_min[j],interval_max[j],num_par[j]))
param_domain=[]
for b in range(num_par[0]):
for d in range(num_par[1]):
param_domain.append([order,grids[0][b],grids[1][d]])
Parallel(n_jobs=num_cores)(delayed(run_single_test)(cand, param_domain[cand]) for cand in range(len(param_domain)))
#! /bin/bash
test=$1
basis=$2
for mesh in 50 100 200 400
do
./../../bin1D/main_dec.out $mesh
mv ./TEST/sol_prim_last.dat ./TEST/test_$test"_"$basis"_Mood_"$mesh".dat"
mv ./TEST/FluxIndicator_last.dat ./TEST/test_$test"_"$basis"_FluxIndicator_"$mesh".dat"
mv ./TEST/DiagIndicator_last.dat ./TEST/test_$test"_"$basis"_DiagIndicator_"$mesh".dat"
rm -rf ./TEST/sol_prim*
rm -rf ./TEST/sol_cons*
rm -rf ./TEST/Flux*
rm -rf ./TEST/Entrop*
rm -rf ./TEST/Diag*
done
#! /bin/bash
test=$1
basis=$2
for mesh in 50 100 200 400
do
./../../bin1D/main_dec.out $mesh
mv ./TEST/sol_prim_last.dat ./TEST/test$test"_"$basis"_Mood_"$mesh".dat"
mv ./TEST/FluxIndicator_last.dat ./TEST/test$test"_"$basis"_FluxIndicator_"$mesh".dat"
mv ./TEST/DiagIndicator_last.dat ./TEST/test$test"_"$basis"_DiagIndicator_"$mesh".dat"
rm -rf ./TEST/sol_prim*
rm -rf ./TEST/sol_cons*
rm -rf ./TEST/Flux*
rm -rf ./TEST/Entrop*
rm -rf ./TEST/Diag*
done
100 nt
1 itype 1: P1, 2: B2, 3: P2, 4:P3, 5: B3
2 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??
1. alpha_jump: parameter used in Burman stabilization term
1. alpha_jump: parameter used in Burman stabilization term
0.1 cfl
100000 ktmax
3.0 tmax
100 ifre
1 test case Euler: 0 - smooth (convergence), 1 - Sod, 2 - Shu-Osher, 3 - Woodward-Colella, 6 - Woodward-Colella left
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
testnr=3
conv=3
orders=["1","2","3"]
styles = ['--^', '-s', ':o']
msize=[7,6,7]
meve=[18,12,6]
lwd=[1,1,2.5]
vars=["rho","v","p" ]#, "Momentum", "Energy"]
#limits=[[[0.45,0.65], [0.37, 0.47]]] Sod first area
#limits=[[[0.63,0.83], [0.1, 0.3]]] #sod second area
limits=[[[0.,1], [0.1, 1.1]], [[0.,1.],[-0.05,1.4]], [[0.,1.],[0.07,1.1]]]
mesh=["200"]
#or N in [int(2**k) for k in range()]:
sols=[]
ex_sol=[]
dashes = [10, 5, 100, 5]
for i, N in enumerate(mesh):
for var, varname in enumerate(vars):
fig, ax = plt.subplots()
#plt.title(varname+", Shu-Osher test, N="+str(N))
#axins = zoomed_inset_axes(ax, 2.5, loc=3)
for ord, ordname in enumerate(orders):
sols.append(np.loadtxt("./B"+ordname+"/TEST/test_Sod_B"+ordname+"_Mood_"+N+".dat"))
ax.plot(sols[ord][:,0],sols[ord][:,var+1],styles[ord],linewidth=lwd[ord],markersize=msize[ord],markevery=meve[ord], label="B"+ordname+", 400 DoFs")
#axins.plot(sols[ord][:,0],sols[ord][:,var+1])
ex_sol.append(np.loadtxt("Reference_Sod/test_Sod_exact_"+varname+".dat"))
ax.plot(ex_sol[var][:,0],ex_sol[var][:,1],'k',label="Exact")
ax.set_xlim(limits[var][0])
ax.set_ylim(limits[var][1])
ax.legend(loc=1)
ax.legend(frameon=False)
#mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
plt.savefig("Test_1D_Sod_"+varname+"_"+N+".pdf")
# plt.show()
plt.clf()
plt.close
limits=[[[0.45,0.85], [0.1, 0.46]], [[0.45,0.85],[-0.05,1.4]], [[0.45,0.85],[0.07,0.43]]] #sod second area
mesh=["200"]
#or N in [int(2**k) for k in range()]:
sols=[]
ex_sol=[]
for i, N in enumerate(mesh):
for var, varname in enumerate(vars):
fig, ax = plt.subplots()
for ord, ordname in enumerate(orders):
sols.append(np.loadtxt("./B"+ordname+"/TEST/test_Sod_B"+ordname+"_Mood_"+N+".dat"))
ax.plot(sols[ord][:,0],sols[ord][:,var+1],styles[ord],linewidth=lwd[ord],markersize=msize[ord],markevery=meve[ord],label="B"+ordname+", 400 DoFs")
ex_sol.append(np.loadtxt("Reference_Sod/test_Sod_exact_"+varname+".dat"))
ax.plot(ex_sol[var][:,0],ex_sol[var][:,1],'k',label="Exact")
ax.set_xlim(limits[var][0])
ax.set_ylim(limits[var][1])
ax.legend(loc=1)
ax.legend(frameon=False)
#mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5")
plt.savefig("Test_1D_Sod_"+varname+"_"+N+"_zoom.pdf")
# plt.show()
plt.clf()
plt.close
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
testnr=3
conv=3
orders=["1","2","3"]
styles = ['--^', '-s', ':o']
msize=[7,6,7]
meve=[18,12,6]