diff --git a/.gitignore b/.gitignore index 7b2ef8c..32d0480 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +<<<<<<< HEAD +# Python files +======= # MatCL files matlab/*.mexa64 matlab/*.mexmaci64 @@ -8,6 +11,7 @@ matlab/*.mexw64 */*.m~ # python files +>>>>>>> 98f26b21b4 __pycache__/ *.py[cod] *$py.class diff --git a/examples/ideal_MHD.py b/examples/ideal_MHD.py new file mode 100644 index 0000000..fec3c65 --- /dev/null +++ b/examples/ideal_MHD.py @@ -0,0 +1,87 @@ +import numpy as np +import sys +import math + + + +prepare_vars() + + +N_fac = 1 + +I_Mesh['NODES_X'] = int(N_fac * 58) +I_Mesh['NODES_Y'] = int(N_fac * 100) +I_Mesh['NODES_Z'] = int(N_fac * 12) + +I_Mesh['XMIN'] = 0.0 +I_Mesh['YMIN'] = 0.0 +I_Mesh['ZMIN'] = 0.0 + +I_Mesh['XMAX'] = 2.0/math.sqrt(3) +I_Mesh['YMAX'] = 2.0 +I_Mesh['ZMAX'] = 1.0 + +I_TI['final_time'] = 5 +I_TI['cfl'] = 0.3 + +dt = I_TI['cfl'] * 2.0 / float(I_Mesh['NODES_Y']) +num_steps = math.ceil(I_TI['final_time']/dt) +dt = I_TI['final_time'] / num_steps + +I_TI['time_integrator'] = 'CarpenterKennedy2N54' + +I_TI['DT'] = dt +I_TI['num_steps'] = num_steps + +# I_Tech['device'] = 1 + +I_Tech['REAL'] = 'double' +I_Tech['REAL4'] = I_Tech['REAL'] + "4" +I_Tech['memory_layout'] = 'USE_ARRAY_OF_STRUCTURES'#'USE_STRUCTURE_OF_ARRAYS' + + +I_BalanceLaws['NUM_CONSERVED_VARS'] = 8 +I_BalanceLaws['NUM_AUXILIARY_VARS'] = 4 +I_BalanceLaws['NUM_TOTAL_VARS'] = I_BalanceLaws['NUM_CONSERVED_VARS'] + I_BalanceLaws['NUM_AUXILIARY_VARS'] + +# Compiler based optimizations + +if I_Tech['REAL'] == 'float': + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant' +else: + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only' + +I_RunOps['periodicx'] = 'None' #'USE_PERIODIC_X' +I_RunOps['periodicy'] = 'None' #'USE_PERIODIC_Y' +I_RunOps['periodicz'] = 'None' #'USE_PERIODIC_Z' + +I_RunOps['order'] = 4 +I_RunOps['conservation_laws'] = 'ideal_MHD' +I_RunOps['testcase'] = 'alfven_periodic' +I_RunOps['plot_numerical_solution'] = 'z' +I_RunOps['save_integrals_over_time'] = False +I_RunOps['norm'] = 'LInf' + +## Initialize + +field_u1, field_u2 = initialize() + + +print('Testcase: ' + I_RunOps['testcase'] + '\nOrder: ' + str(I_RunOps['order']) + '\nTime integrator: ' + \ + I_TI['time_integrator'] +' \nDT: ' + str(I_TI['DT']) + '\nN_STEPS: ' +str(I_TI['num_steps']) + " FINAL_TIME: " + \ + str(I_TI['final_time']) + '\nDX: ' + str(I_Mesh['DX']) + ' NODES_X: ' +str(I_Mesh['NODES_X'])+'\nDY: '+ str(I_Mesh['DY']) + \ + ' NODES_Y: ' + str(I_Mesh['NODES_Y']) + '\nDZ: ' + str(I_Mesh['DZ'])+ ' NODES_Z: ' + str(I_Mesh['NODES_Z']) + '\nREAL: ' + \ + str(I_Tech['REAL'])) + +compute_numerical_solution(field_u1, field_u2) + +rel_err = I_Results['rel_err'] +for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + print("Relative Error of Field Component " + str(comp) +": " + str(100*rel_err[comp])) + + +field_u1_reshaped = np.reshape(field_u1, (I_Tech['NUM_NODES_PAD'], I_BalanceLaws['NUM_TOTAL_VARS'])) + + +save_all_variables(field_u1_reshaped, "output.hdf5", ['rho', 'px', 'py', 'pz', 'E', 'Bx', 'By', 'Bz']) +plot_2D(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', 1,1) diff --git a/examples/ideal_MHD_B0.py b/examples/ideal_MHD_B0.py new file mode 100644 index 0000000..62a6b9d --- /dev/null +++ b/examples/ideal_MHD_B0.py @@ -0,0 +1,90 @@ +import numpy as np +import sys +import math + +prepare_vars() + + +N_fac = 1 + +I_Mesh['NODES_X'] = int(N_fac * 100) +I_Mesh['NODES_Y'] = int(N_fac * 100) +I_Mesh['NODES_Z'] = int(N_fac * 100) + +I_Mesh['XMIN'] = 0.0 +I_Mesh['YMIN'] = 0.0 +I_Mesh['ZMIN'] = 0.0 + +I_Mesh['XMAX'] = 1.0 +I_Mesh['YMAX'] = 1.0 +I_Mesh['ZMAX'] = 1.0 + +I_TI['final_time'] = 0.1 +I_TI['cfl'] = 0.2 +k = 100 + +dt = I_TI['cfl'] * 2.0 / float(I_Mesh['NODES_Y']) # this has to be estimated +num_steps = math.ceil(I_TI['final_time']/dt) +dt = I_TI['final_time'] / num_steps + +I_TI['time_integrator'] = 'CarpenterKennedy2N54' + +I_TI['DT'] = dt +I_TI['num_steps'] = num_steps + +# I_Tech['device'] = 1 + +I_Tech['REAL'] = 'double' +I_Tech['REAL4'] = I_Tech['REAL'] + "4" +I_Tech['memory_layout'] = 'USE_ARRAY_OF_STRUCTURES' #'USE_STRUCTURE_OF_ARRAYS' + + +I_BalanceLaws['NUM_CONSERVED_VARS'] = 8 +I_BalanceLaws['NUM_AUXILIARY_VARS'] = 4 +I_BalanceLaws['NUM_TOTAL_VARS'] = I_BalanceLaws['NUM_CONSERVED_VARS'] + I_BalanceLaws['NUM_AUXILIARY_VARS'] + +# Compiler based optimizations + +if I_Tech['REAL'] == 'float': + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant' +else: + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only' + +I_RunOps['periodic'] = 'none'# 'USE_PERIODIC' +I_RunOps['vr'] = 'none'#'USE_VR_KUSANO' + +I_RunOps['order'] = 4 +I_RunOps['conservation_laws'] = 'ideal_MHD' +I_RunOps['testcase'] = 'B0' +I_RunOps['plot_numerical_solution'] = 'z' +I_RunOps['save_integrals_over_time'] = False +I_RunOps['norm'] = 'LInf' + +## Initialize + +field_u1, field_u2 = initialize() + + +print('Testcase: ' + I_RunOps['testcase'] + '\nOrder: ' + str(I_RunOps['order']) + '\nTime integrator: ' + \ + I_TI['time_integrator'] +' \nDT: ' + str(I_TI['DT']) + '\nN_STEPS: ' +str(I_TI['num_steps']) + " FINAL_TIME: " + \ + str(I_TI['final_time']) + '\nDX: ' + str(I_Mesh['DX']) + ' NODES_X: ' +str(I_Mesh['NODES_X'])+'\nDY: '+ str(I_Mesh['DY']) + \ + ' NODES_Y: ' + str(I_Mesh['NODES_Y']) + '\nDZ: ' + str(I_Mesh['DZ'])+ ' NODES_Z: ' + str(I_Mesh['NODES_Z']) + '\nREAL: ' + \ + str(I_Tech['REAL'])) + +print(str(k) + " snapshots") +for i in range(k): + compute_numerical_solution(field_u1, field_u2) + save_all_variables(field_u1, "results/output" + str(i) + ".hdf5", ['rho', 'px', 'py', 'pz', 'E', 'Bx', 'By', 'Bz']) + print(i) + +#rel_err = I_Results['rel_err'] +#for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + # print("Relative Error of Field Component " + str(comp) +": " + str(100*rel_err[comp])) + + +#field_u1_reshaped = np.reshape(field_u1, (I_Tech['NUM_NODES_PAD'], I_BalanceLaws['NUM_TOTAL_VARS'])) + + + +#plot_fieldlines_B(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', Field_Bx, Field_By) +#plot_2D(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', Field_By,Field_By) diff --git a/examples/ideal_MHD_dipole.py b/examples/ideal_MHD_dipole.py new file mode 100644 index 0000000..b11c8aa --- /dev/null +++ b/examples/ideal_MHD_dipole.py @@ -0,0 +1,98 @@ +import numpy as np +import sys +import math + +prepare_vars() + + +N_fac = 1 + +I_Mesh['NODES_X'] = int(N_fac * 100) +I_Mesh['NODES_Y'] = int(N_fac * 100) +I_Mesh['NODES_Z'] = int(N_fac * 100) + +I_Mesh['XMIN'] = 0.0 +I_Mesh['YMIN'] = 0.0 +I_Mesh['ZMIN'] = 0.0 + +I_Mesh['XMAX'] = 1.0 +I_Mesh['YMAX'] = 1.0 +I_Mesh['ZMAX'] = 1.0 + +I_TI['final_time'] = 0.1 +I_TI['cfl'] = 0.2 +k = 10 + +dt = I_TI['cfl'] * 2.0 / float(I_Mesh['NODES_Y']) # this has to be estimated +num_steps = math.ceil(I_TI['final_time']/dt) +dt = I_TI['final_time'] / num_steps + +I_TI['time_integrator'] = 'CarpenterKennedy2N54' + +I_TI['DT'] = dt +I_TI['num_steps'] = num_steps + +# I_Tech['device'] = 1 + +I_Tech['REAL'] = 'float' #'double' +I_Tech['REAL4'] = I_Tech['REAL'] + "4" +I_Tech['memory_layout'] = 'USE_ARRAY_OF_STRUCTURES' #'USE_STRUCTURE_OF_ARRAYS' + + +I_BalanceLaws['NUM_CONSERVED_VARS'] = 8 +I_BalanceLaws['NUM_AUXILIARY_VARS'] = 4 +I_BalanceLaws['NUM_TOTAL_VARS'] = I_BalanceLaws['NUM_CONSERVED_VARS'] + I_BalanceLaws['NUM_AUXILIARY_VARS'] + +# Compiler based optimizations + +if I_Tech['REAL'] == 'float': + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant' +else: + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only' + +I_RunOps['periodicx'] = 'none'# 'USE_PERIODIC' +I_RunOps['periodicy'] = 'none' +I_RunOps['periodicz'] = 'none' + +I_RunOps['vr'] = 'none'#'USE_VR_KUSANO' + +I_RunOps['order'] = 2 +I_RunOps['conservation_laws'] = 'ideal_MHD' +I_RunOps['testcase'] = 'far_dipole' +I_RunOps['plot_numerical_solution'] = 'z' +I_RunOps['save_integrals_over_time'] = False +I_RunOps['norm'] = 'LInf' + +## Initialize + +field_u1, field_u2 = initialize() + +if len(sys.argv) > 1: + field_u1 = reload_all_variables(field_u1, sys.argv[1], ['rho', 'px', 'py', 'pz', 'E', 'Bx', 'By', 'Bz']) + print(f"reloaded the conserved variables out of {sys.argv[1]}") + +print('Testcase: ' + I_RunOps['testcase'] + '\nOrder: ' + str(I_RunOps['order']) + '\nTime integrator: ' + \ + I_TI['time_integrator'] +' \nDT: ' + str(I_TI['DT']) + '\nN_STEPS: ' +str(I_TI['num_steps']) + " FINAL_TIME: " + \ + str(I_TI['final_time']) + '\nDX: ' + str(I_Mesh['DX']) + ' NODES_X: ' +str(I_Mesh['NODES_X'])+'\nDY: '+ str(I_Mesh['DY']) + \ + ' NODES_Y: ' + str(I_Mesh['NODES_Y']) + '\nDZ: ' + str(I_Mesh['DZ'])+ ' NODES_Z: ' + str(I_Mesh['NODES_Z']) + '\nREAL: ' + \ + str(I_Tech['REAL'])) + +print(str(k) + " snapshots") +ct = False + +for i in range(k): + ct = compute_numerical_solution(field_u1, field_u2, ct) + save_all_variables(field_u1, "results/output" + str(i) + ".hdf5", ['rho', 'px', 'py', 'pz', 'E', 'Bx', 'By', 'Bz']) + print(i) + +#rel_err = I_Results['rel_err'] +#for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + # print("Relative Error of Field Component " + str(comp) +": " + str(100*rel_err[comp])) + + +#field_u1_reshaped = np.reshape(field_u1, (I_Tech['NUM_NODES_PAD'], I_BalanceLaws['NUM_TOTAL_VARS'])) + + + +#plot_fieldlines_B(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', Field_Bx, Field_By) +#plot_2D(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', Field_By,Field_By) diff --git a/examples/ideal_MHD_ulf_p.py b/examples/ideal_MHD_ulf_p.py new file mode 100644 index 0000000..eaaad97 --- /dev/null +++ b/examples/ideal_MHD_ulf_p.py @@ -0,0 +1,98 @@ +import numpy as np +import sys +import math + +prepare_vars() + + +N_fac = 1 + +I_Mesh['NODES_X'] = int(N_fac * 100) +I_Mesh['NODES_Y'] = int(N_fac * 100) +I_Mesh['NODES_Z'] = int(N_fac * 100) + +I_Mesh['XMIN'] = 0.0 +I_Mesh['YMIN'] = 0.0 +I_Mesh['ZMIN'] = 0.0 + +I_Mesh['XMAX'] = 1.0 +I_Mesh['YMAX'] = 1.0 +I_Mesh['ZMAX'] = 1.0 + +I_TI['final_time'] = 0.1 +I_TI['cfl'] = 0.2 +k = 10 + +dt = I_TI['cfl'] * 2.0 / float(I_Mesh['NODES_Y']) # this has to be estimated +num_steps = math.ceil(I_TI['final_time']/dt) +dt = I_TI['final_time'] / num_steps + +I_TI['time_integrator'] = 'CarpenterKennedy2N54' + +I_TI['DT'] = dt +I_TI['num_steps'] = num_steps + +# I_Tech['device'] = 1 + +I_Tech['REAL'] = 'float' #'double' +I_Tech['REAL4'] = I_Tech['REAL'] + "4" +I_Tech['memory_layout'] = 'USE_ARRAY_OF_STRUCTURES' #'USE_STRUCTURE_OF_ARRAYS' + + +I_BalanceLaws['NUM_CONSERVED_VARS'] = 8 +I_BalanceLaws['NUM_AUXILIARY_VARS'] = 4 +I_BalanceLaws['NUM_TOTAL_VARS'] = I_BalanceLaws['NUM_CONSERVED_VARS'] + I_BalanceLaws['NUM_AUXILIARY_VARS'] + +# Compiler based optimizations + +if I_Tech['REAL'] == 'float': + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant' +else: + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only' + +I_RunOps['periodicx'] = 'none'# 'USE_PERIODIC' +I_RunOps['periodicy'] = 'none' +I_RunOps['periodicz'] = 'none' + +I_RunOps['vr'] = 'none'#'USE_VR_KUSANO' + +I_RunOps['order'] = 6 +I_RunOps['conservation_laws'] = 'ideal_MHD' +I_RunOps['testcase'] = 'ufp' +I_RunOps['plot_numerical_solution'] = 'z' +I_RunOps['save_integrals_over_time'] = False +I_RunOps['norm'] = 'LInf' + +## Initialize + +field_u1, field_u2 = initialize() + +if len(sys.argv) > 1: + field_u1 = reload_all_variables(field_u1, sys.argv[1], ['rho', 'px', 'py', 'pz', 'E', 'Bx', 'By', 'Bz']) + print(f"reloaded the conserved variables out of {sys.argv[1]}") + +print('Testcase: ' + I_RunOps['testcase'] + '\nOrder: ' + str(I_RunOps['order']) + '\nTime integrator: ' + \ + I_TI['time_integrator'] +' \nDT: ' + str(I_TI['DT']) + '\nN_STEPS: ' +str(I_TI['num_steps']) + " FINAL_TIME: " + \ + str(I_TI['final_time']) + '\nDX: ' + str(I_Mesh['DX']) + ' NODES_X: ' +str(I_Mesh['NODES_X'])+'\nDY: '+ str(I_Mesh['DY']) + \ + ' NODES_Y: ' + str(I_Mesh['NODES_Y']) + '\nDZ: ' + str(I_Mesh['DZ'])+ ' NODES_Z: ' + str(I_Mesh['NODES_Z']) + '\nREAL: ' + \ + str(I_Tech['REAL'])) + +print(str(k) + " snapshots") +ct = False + +for i in range(k): + ct = compute_numerical_solution(field_u1, field_u2, ct) + save_all_variables(field_u1, "results/output" + str(i) + ".hdf5", ['rho', 'px', 'py', 'pz', 'E', 'Bx', 'By', 'Bz']) + print(i) + +#rel_err = I_Results['rel_err'] +#for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + # print("Relative Error of Field Component " + str(comp) +": " + str(100*rel_err[comp])) + + +#field_u1_reshaped = np.reshape(field_u1, (I_Tech['NUM_NODES_PAD'], I_BalanceLaws['NUM_TOTAL_VARS'])) + + + +#plot_fieldlines_B(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', Field_Bx, Field_By) +#plot_2D(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', Field_By,Field_By) diff --git a/examples/linear_constant_advection.py b/examples/linear_constant_advection.py new file mode 100644 index 0000000..e210597 --- /dev/null +++ b/examples/linear_constant_advection.py @@ -0,0 +1,80 @@ +import numpy as np +import sys +import math + + + +prepare_vars() + + +N = 500 +I_Mesh['NODES_X'] = N +I_Mesh['NODES_Y'] = N +I_Mesh['NODES_Z'] = N + +I_Mesh['XMIN'] = -1.0 +I_Mesh['YMIN'] = -1.0 +I_Mesh['ZMIN'] = -1.0 + +I_Mesh['XMAX'] = 1.0 +I_Mesh['YMAX'] = 1.0 +I_Mesh['ZMAX'] = 1.0 + +I_TI['final_time'] = 1 +I_TI['cfl'] = 0.95 + +dt = I_TI['cfl'] * 2.0 / float(N) +num_steps = math.ceil(I_TI['final_time']/dt) +dt = I_TI['final_time'] / num_steps + +I_TI['time_integrator'] = 'CarpenterKennedy2N54' + +I_TI['DT'] = dt +I_TI['num_steps'] = num_steps + +# I_Tech['device'] = 1 + +I_Tech['REAL'] = 'double' +I_Tech['REAL4'] = I_Tech['REAL'] + '4' + +I_BalanceLaws['NUM_CONSERVED_VARS'] = 1 +I_BalanceLaws['NUM_AUXILILIARY_VARS'] = 0 +I_BalanceLaws['NUM_TOTAL_VARS'] = I_BalanceLaws['NUM_CONSERVED_VARS'] + I_BalanceLaws['NUM_AUXILIARY_VARS'] + +# Compiler based optimizations + +if I_Tech['REAL'] == 'float': + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant' +else: + I_Tech['optimizations'] = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only' + +I_RunOps['periodic'] = 'USE_PERIODIC' + +I_RunOps['order'] = 4 +I_RunOps['conservation_laws'] = 'linear_constant_advection' +I_RunOps['testcase'] = 'donut' +I_RunOps['plot_numerical_solution'] = 'z' +I_RunOps['save_integrals_over_time'] = False +I_RunOps['norm'] = 'LInf' + +## Initialize + +field_u1, field_u2 = initialize() + + +print('Testcase: ' + I_RunOps['testcase'] + '\nOrder: ' + str(I_RunOps['order']) + '\nTime integrator: ' + \ + I_TI['time_integrator'] +' \nDT: ' + str(I_TI['DT']) + '\nN_STEPS: ' +str(I_TI['num_steps']) + " FINAL_TIME: " + \ + str(I_TI['final_time']) + '\nDX: ' + str(I_Mesh['DX']) + ' NODES_X: ' +str(I_Mesh['NODES_X'])+'\nDY: '+ str(I_Mesh['DY']) + \ + ' NODES_Y: ' + str(I_Mesh['NODES_Y']) + '\nDZ: ' + str(I_Mesh['DZ'])+ ' NODES_Z: ' + str(I_Mesh['NODES_Z']) + '\nREAL: ' + \ + str(I_Tech['REAL'])) + +compute_numerical_solution(field_u1, field_u2) + +rel_err = I_Results['rel_err'] +for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + print("Relative Error of Field Component " + str(comp) +": " + str(100*rel_err[comp])) + + +field_u1_reshaped = np.reshape(field_u1, (I_Tech['num_nodes_pad'], I_BalanceLaws['NUM_CONSERVED_VARS'])) + +plot_2D(field_u1_reshaped, I_RunOps['plot_numerical_solution'], I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z'], 'Numerical Solution', 1, 1) diff --git a/examples/starter.sh b/examples/starter.sh new file mode 100755 index 0000000..5b6cea4 --- /dev/null +++ b/examples/starter.sh @@ -0,0 +1,10 @@ +#!/bin/sh + +cat ../python/*.py > runfile.py +cat ../python/BalanceLaws/*.py >> runfile.py +cp ../python/cl/cl.py cl.py +cat $1 >> runfile.py +python3 runfile.py $2 +rm cl.py +#rm runfile.py + diff --git a/include_physics/ideal_MHD.h b/include_physics/ideal_MHD.h index f4ca570..8cc07c1 100644 --- a/include_physics/ideal_MHD.h +++ b/include_physics/ideal_MHD.h @@ -406,6 +406,7 @@ inline void calc_flux_f(const REAL *u, REAL *flux) { return; } + inline void calc_flux_g(const REAL *u, REAL *flux){ REAL BB = u[Field_Bx]*u[Field_Bx] + u[Field_By]*u[Field_By] + u[Field_Bz]*u[Field_Bz]; flux[Field_rho] = u[Field_rho_uy]; @@ -434,6 +435,7 @@ inline void calc_flux_h(const REAL *u, REAL *flux){ + #if defined USE_BOUNDARY_FLUX_HLL inline void calc_num_flux(REAL al, REAL ar, REAL *ul, REAL *ur, REAL *fluxl, REAL *fluxr, REAL *num_flux){ @@ -453,6 +455,7 @@ inline void calc_num_flux(REAL al, REAL ar, REAL *ul, REAL *ur, REAL *fluxl, REA } + inline void calc_hll_speeds(REAL* ul, REAL* ur, REAL *cl, REAL* cr, int dir) { #ifdef USE_HLL_SPEED_BKW // Speeds from @@ -609,6 +612,7 @@ inline void add_surface_terms(REAL time, uint ix, uint iy, uint iz, const global #ifndef USE_PERIODIC_X if (check_bound_xr(ix, 1)) { calc_hll_speeds(um, ub, &alx, &arx, 0); + calc_flux_f(um, fluxm); calc_flux_f(ub, fluxb); calc_num_flux(alx, arx, um, ub, fluxm, fluxb, flux); diff --git a/include_testcases/ideal_MHD_far_dipole.h b/include_testcases/ideal_MHD_far_dipole.h index eec1652..f3c211f 100644 --- a/include_testcases/ideal_MHD_far_dipole.h +++ b/include_testcases/ideal_MHD_far_dipole.h @@ -9,7 +9,7 @@ // A dipole outside of the domain // paramters REAL4 constant M = {0.0, 0.0, 1.0, 0.0}; -REAL4 constant X_0 = {-1.0, 0.5, 0.0, 0.0}; +REAL4 constant X_0 = {-1.0, 0.5, 0.5, 0.0}; REAL constant A0 = 0.1; REAL constant omega = 6.28 * 10; /* diff --git a/python/BalanceLaws/compute_numerical_solution.py b/python/BalanceLaws/compute_numerical_solution.py new file mode 100644 index 0000000..ae37055 --- /dev/null +++ b/python/BalanceLaws/compute_numerical_solution.py @@ -0,0 +1,143 @@ +def compute_numerical_solution(field_u1, field_u2, c_time = 'none'): + global I_mesh, I_TI, I_BalanceLaws, I_Tech, I_RunOps, I_Results + + + if I_Tech['REAL'] == 'float': + current_time = np.zeros(2, dtype=np.float32) + components = np.zeros(2, dtype=np.float32) + norm_output = np.zeros((I_Tech['num_groups'], 1), dtype=np.float32) + Lerror = np.zeros((I_TI['num_steps']+1, I_BalanceLaws['NUM_CONSERVED_VARS']), dtype=np.float32) + else: + current_time = np.zeros(2, dtype=np.float64) + components = np.zeros(2, dtype=np.float64) + norm_output = np.zeros((I_Tech['num_groups'], 1), dtype=np.float64) + Lerror = np.zeros((I_TI['num_steps']+1, I_BalanceLaws['NUM_CONSERVED_VARS']), dtype=np.float64) + + if type(c_time) == type(np.array([0.0, 0.0])): + current_time = c_time + + if I_TI['time_integrator'] == 'CarpenterKennedy2N54': + RK_Step_a = ['CarpenterKennedy2N54_1a', 'CarpenterKennedy2N54_2a', \ + 'CarpenterKennedy2N54_3a', 'CarpenterKennedy2N54_4a', \ + 'CarpenterKennedy2N54_5a'] + + RK_Step_b = ['CarpenterKennedy2N54_1b', 'CarpenterKennedy2N54_2b', \ + 'CarpenterKennedy2N54_3b', 'CarpenterKennedy2N54_4b', \ + 'CarpenterKennedy2N54_5b'] + RK_Step = [] + for i in range(len(RK_Step_a)): + RK_Step.append('calc_auxiliary_vars_1_2_args') + RK_Step.append(RK_Step_a[i]) + RK_Step.append(RK_Step_b[i]) + + time_integrator_num_fields = 2; + + if I_TI['time_integrator'] == 'SSPRK33': + RK_Step = [] + RK_Step.append('calc_auxiliary_vars_1_3_args') + RK_Step.append('SSPRK33_1') + RK_Step.append('SSPRK33_2a') + RK_Step.append('calc_auxiliary_vars_2_3_args') + RK_Step.append('SSPRK33_2b') + RK_Step.append('SSPRK33_3a') + RK_Step.append('calc_auxiliary_vars_3_3_args') + RK_Step.append('SSPRK33_3b') + RK_Step.append('calc_time_3_args') + time_integrator_num_fields = 3 + + num_steps_run = I_TI['num_steps'] + kernel_runtime = 0 + RK_block_size = 15000 + + if time_integrator_num_fields == 2: + if I_RunOps['save_integrals_over_time']: + for step in range(I_TI['num_steps']): + for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + components[0] = comp + cl.run_kernel('analytical_u', I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u2, current_time) + norm_output[:] = 0 + if I_RunOps['norm'] == 'L2': + cl.run_kernel('norm2_diff'. I_Tech['g_range'], I_Tech['l_range'], field_u1, field_u2, norm_output, components) + Lerror[comp, step] = sqrt(np.sum(norm_output)) + elif I_RunOps['norm'] == 'LInf': + cl.run_kernel('norm_infty_diff', I_Tech['g_range'], I_Tech['l_range'], field_u1, field_u2, norm_output, components) + Lerror[comp, step] = np.max(norm_output) + cl.run_kernel(RK_Step, I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1, field_u2, current_time) + else: + while num_steps_run > RK_block_size: + kernel_list = RK_Step * RK_block_size + t = cl.run_kernel(kernel_list, I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1, field_u2, current_time) + kernel_runtime = kernel_runtime + t + num_steps_run = num_steps_run - RK_block_size + if num_steps_run > 0: + kernel_list = num_steps_run* RK_Step + t = cl.run_kernel(kernel_list, I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1, field_u2, current_time) + kernel_runtime = kernel_runtime + t + + + if time_integrator_num_fields == 3: + num_nodes = I_Mesh['NODES_X']*I_Mesh['NODES_Y']*I_Mesh['NODES_Z'] + if I_Tech['REAL'] == 'float': + field_u3 = np.zeros((num_nodes*I_BalanceLaws['NUM_TOTAL_VARS'], 1), dtype = np.float32) + else: + field_u3 = np.zeros((num_nodes*I_BalanceLaws['NUM_TOTAL_VARS'], 1), dtype = np.float64) + + if I_RunOps['save_integrals_over_time']: + for step in range(TI['num_steps']): + for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + components[0] = comp + cl.run_kernel('analytical_u', I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u2, current_time) + norm_output[:] = 0 + if I_RunOps['norm'] == 'L2': + cl.run_kernel('norm2_diff'. I_Tech['g_range'], I_Tech['l_range'], field_u1, field_u2, field_u3,norm_output, components) + Lerror[comp, step] = sqrt(np.sum(norm_output)) + elif I_RunOps['norm'] == 'LInf': + cl.run_kernel('norm_infty_diff', I_Tech['g_range'], I_Tech['l_range'], field_u1, field_u2, field_u3,norm_output, components) + Lerror[comp, step] = np.max(norm_output) + cl.run_kernel(RK_Step, I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1, field_u2, field_u3, current_time) + + else: + while num_steps_run > RK_block_size: + kernel_list = RK_Step * RK_block_size + t = cl.run_kernel(kernel_list, I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1, field_u2, field_u3, current_time) + kernel_runtime = kernel_runtime + t + num_steps_run = num_steps_run - RK_block_size + if num_steps_run > 0: + kernel_list = num_steps_run* RK_Step + t = cl.run_kernel(kernel_list, I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1, field_u2, field_u3, current_time) + kernel_runtime = kernel_runtime + t + + + if I_RunOps['save_fields']: + I_Results['field_u'] = field_u1 + + print(f"current_time[0] = {current_time[0]} current_time[1] = {current_time[1]}") + current_time[1] = I_TI['final_time'] + cl.run_kernel('analytical_u', I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u2, current_time) + abs_err = np.zeros((I_BalanceLaws['NUM_CONSERVED_VARS'],1), dtype = np.float64) + rel_err = np.zeros(( I_BalanceLaws['NUM_CONSERVED_VARS'],1),dtype=np.float64) + + for comp in range(I_BalanceLaws['NUM_CONSERVED_VARS']): + components[0] = comp + norm_output[:]= 0 + if I_RunOps['norm'] == 'L2': +# cl.run_kernel('norm2_diff', I_Tech['g_range'], I_Tech['l_range'],field_u1, field_u2, norm_output, components) + abs_err[comp] = math.sqrt(np.sum(norm_output)) + norm_output[:] = 0 + # cl.run_kernel('norm2', I_Tech['g_range'], I_Tech['l_range'], field_u2, norm_output, components) + rel_err[comp] = abs_err[comp] / math.sqrt(np.sum(norm_output)) + elif I_RunOps['norm'] == 'LInf': +# cl.run_kernel('norm_infty_diff', I_Tech['g_range'], I_Tech['l_range'],field_u1, field_u2, norm_output, components) + abs_err[comp] = np.max(norm_output) + norm_output[:] = 0 + # cl.run_kernel('norm_infty', I_Tech['g_range'], I_Tech['l_range'], field_u2, norm_output, components) + rel_err[comp] = abs_err[comp] / np.max(norm_output) + + I_Results['abs_err'] = abs_err + I_Results['rel_err'] = rel_err + + if I_RunOps['save_integrals_over_time']: + Lerror[I_TI['num_steps'],:] = I_Results['abs_err'] + I_Results['Lerror_over_time'] = Lerror + I_Results['time'] = np.linspace(0, I_TI['final_time'], len(Lerror)*2) + return current_time diff --git a/python/BalanceLaws/generate_setting.py b/python/BalanceLaws/generate_setting.py new file mode 100644 index 0000000..0f5af59 --- /dev/null +++ b/python/BalanceLaws/generate_setting.py @@ -0,0 +1,23 @@ +def generate_settings(map, keys): + # based on the matlab counterpart + # Atomatically generates a formatted settings string that contains all + # compiler optimizations and OpenCL defines of map with the given keys where + # keys is a list, e.g. settings_mesh = generate_settings(I_Mesh, ['DX', 'DY', 'DZ']) + + settings = ''; + for key in keys: + if key in map: + if not(type(map[key]) is float or type(map[key]) is int): + if key== 'optimizations': + settings = settings + map[key] + elif 'USE' in map[key]: + settings = settings + ' -D' + map[key] + '=1' + else: + settings = settings + ' -D' + key + '=' + str(map[key]) + else: + settings = settings + ' -D' + key + "=" + str(map[key]) + else: + print('Unknown identifier ' + key) + return settings + + diff --git a/python/BalanceLaws/initialize.py b/python/BalanceLaws/initialize.py new file mode 100644 index 0000000..90d7c24 --- /dev/null +++ b/python/BalanceLaws/initialize.py @@ -0,0 +1,157 @@ +import cl + +def initialize(): + global I_Mesh, I_TI, I_BalanceLaws, I_Tech, I_RunOps + + cl.create_ctx() + + # a python float is a 64bit float => float means double + if I_RunOps['periodicx'] == 'USE_PERIODIC_X': + DX = float(I_Mesh['XMAX'] - I_Mesh['XMIN']) / float(I_Mesh['NODES_X']) + else: + DX = float(I_Mesh['XMAX'] - I_Mesh['XMIN']) / (float(I_Mesh['NODES_X'])-1) + + if I_RunOps['periodicy'] == 'USE_PERIODIC_Y': + DY = float(I_Mesh['YMAX'] - I_Mesh['YMIN']) / float(I_Mesh['NODES_Y']) + else: + DY = float(I_Mesh['YMAX'] - I_Mesh['YMIN']) / (float(I_Mesh['NODES_Y'])-1) + + if I_RunOps['periodicz'] == 'USE_PERIODIC_Z': + DZ = float(I_Mesh['ZMAX'] - I_Mesh['ZMIN']) / float(I_Mesh['NODES_Z']) + else: + DZ = float(I_Mesh['ZMAX'] - I_Mesh['ZMIN']) / (float(I_Mesh['NODES_Z'])-1) + + + I_Mesh['DX'] = DX; I_Mesh['DY'] = DY; I_Mesh['DZ'] = DZ + + I_BalanceLaws['g_range'] = np.array([I_Mesh['NODES_X'], I_Mesh['NODES_Y'], I_Mesh['NODES_Z']], dtype=np.uint32) + # I_BalanceLaws['g_range'] = np.array([I_Mesh['NODES_Z'], I_Mesh['NODES_Y'], I_Mesh['NODES_X']], dtype=np.uint32) + # I_BalanceLaws['g_range'] = np.array([I_Mesh['NODES_X']*I_BalanceLaws['NUM_TOTAL_VARS'],I_Mesh['NODES_Y'],I_Mesh['NODES_Z']], dtype=np.uint32) + + I_BalanceLaws['l_range'] = np.array([0], dtype=np.uint32) + + num_nodes = I_Mesh['NODES_X']*I_Mesh['NODES_Y']*I_Mesh['NODES_Z'] + + if cl.iscpu: + group_size=16 + else: + if (num_nodes > cl.group_size): + group_size = cl.group_size + else: + group_size = 2**floor(log(num_nodes) / log(2)) + + + num_nodes_pad = math.ceil(num_nodes/group_size)*group_size + num_groups = math.ceil(num_nodes_pad/group_size) + + I_Tech['NUM_NODES_PAD'] = int(num_nodes_pad) + I_Tech['num_groups'] = int(num_groups) + I_Tech['W_SIZE'] = int(group_size) + + I_Tech['g_range'] = [I_Tech['NUM_NODES_PAD'],1,1] + I_Tech['l_range' ] = [I_Tech['W_SIZE'],1,1] + + + if I_Tech['REAL'] == 'float': + field_u1 = np.zeros((num_nodes_pad*I_BalanceLaws['NUM_TOTAL_VARS'], 1), dtype=np.float32) + field_u2 = np.zeros((num_nodes_pad*I_BalanceLaws['NUM_TOTAL_VARS'], 1), dtype=np.float32) + else: + field_u1 = np.zeros((num_nodes_pad*I_BalanceLaws['NUM_TOTAL_VARS'], 1), dtype = np.float64) + field_u2 = np.zeros((num_nodes_pad*I_BalanceLaws['NUM_TOTAL_VARS'], 1), dtype = np.float64) + + + x = np.linspace(I_Mesh['XMIN'], I_Mesh['XMAX'], I_Mesh['NODES_X']); + y = np.linspace(I_Mesh['YMIN'], I_Mesh['YMAX'], I_Mesh['NODES_Y']); + z = np.linspace(I_Mesh['ZMIN'], I_Mesh['ZMAX'], I_Mesh['NODES_Z']); + + kernel_path_list = [] + + if I_RunOps['order'] == 2: + if I_RunOps['operator_form'] == 'extended': + kernel_path_list.append('../include/2ndOrderExtended.h') + else: # classical operators + kernel_path_list.append('../include/2ndOrder.h') + + elif I_RunOps['order'] == 4: + if I_RunOps['operator_form'] == 'extended': + kernel_path_list.append('../include/4thOrderExtended.h') + else: # classical operators + kernel_path_list.append('../include/4thOrder.h') + + elif I_RunOps['order'] == 6: + if I_RunOps['operator_form'] == 'extended': + kernel_path_list.append('../include/6thOrderExtended.h') + else: # classical operators + kernel_path_list.append('../include/6thOrder.h') + + else: + print('Specify order \n') + + kernel_path_list.append('../include_defines/' + I_RunOps['conservation_laws'] + '_defines.h') + kernel_path_list.append('../include/utils.h') + kernel_path_list.append('../kernel/SBP_operator.cl') + kernel_path_list.append('../include_testcases/' + I_RunOps['conservation_laws'] + '_' + I_RunOps['testcase']+ '.h') + kernel_path_list.append('../include_physics/' + I_RunOps['conservation_laws'] + '.h') + kernel_path_list.append('../kernel/kernel_init.cl') + + settings_tech = generate_settings(I_Tech, ['REAL', 'REAL4', 'optimizations', 'memory_layout', 'NUM_NODES_PAD']) + + settings_mesh = generate_settings(I_Mesh, ['DX', 'DY', 'DZ','NODES_X', 'NODES_Y', 'NODES_Z', 'XMIN', 'XMAX', 'YMIN', 'YMAX','ZMIN', 'ZMAX']) + settings_runops = generate_settings(I_RunOps, ['periodicx', 'periodicy', 'periodicz']) + settings_TI = generate_settings(I_TI, ['DT']) + + settings = settings_tech + settings_mesh + settings_runops + settings_TI + + cl.compile_kernels(kernel_path_list, settings) + + cl.run_kernel( 'init', I_BalanceLaws['g_range'], I_BalanceLaws['l_range'], field_u1) + + + kernel_path_list = [] + + if I_RunOps['order'] == 2: + if I_RunOps['operator_form'] == 'extended': + kernel_path_list.append('../include/2ndOrderExtended.h') + else: # classical operators + kernel_path_list.append('../include/2ndOrder.h') + + elif I_RunOps['order'] == 4: + if I_RunOps['operator_form'] == 'extended': + kernel_path_list.append('../include/4thOrderExtended.h') + else: # classical operators + kernel_path_list.append('../include/4thOrder.h') + + elif I_RunOps['order'] == 6: + if I_RunOps['operator_form'] == 'extended': + kernel_path_list.append('../include/6thOrderExtended.h') + else: # classical operators + kernel_path_list.append('../include/6thOrder.h') + + else: + print('Specify order \n') + + kernel_path_list.append('../include_defines/' + I_RunOps['conservation_laws'] + '_defines.h') + kernel_path_list.append('../include/utils.h') + kernel_path_list.append('../kernel/SBP_operator.cl') + kernel_path_list.append('../include_testcases/' + I_RunOps['conservation_laws'] + '_' + I_RunOps['testcase']+ '.h') + kernel_path_list.append('../include_physics/' + I_RunOps['conservation_laws'] + '.h') + kernel_path_list.append('../kernel/artificial_dissipation.cl') + kernel_path_list.append('../kernel/kernel_time_integrator.cl') + kernel_path_list.append('../kernel/kernel_norm.cl') + + + settings_tech = generate_settings(I_Tech, ['REAL', 'REAL4', 'W_SIZE', 'optimizations', 'memory_layout', 'NUM_NODES_PAD']) + + settings_mesh = generate_settings(I_Mesh, ['DX', 'DY', 'DZ','NODES_X', 'NODES_Y', 'NODES_Z', 'XMIN', 'XMAX', 'YMIN', 'YMAX','ZMIN', 'ZMAX']) + settings_time_integration = generate_settings(I_TI, ['DT']) + settings_runops = generate_settings(I_RunOps, ['periodicx', 'periodicy', 'periodicz']) + settings_TI = generate_settings(I_TI, ['DT']) + + settings = settings_tech + settings_mesh + settings_time_integration +settings_runops + settings_TI + + cl.compile_kernels(kernel_path_list, settings) + + return field_u1, field_u2 + + + diff --git a/python/cl/cl.py b/python/cl/cl.py new file mode 100644 index 0000000..a7d589c --- /dev/null +++ b/python/cl/cl.py @@ -0,0 +1,78 @@ +import pyopencl as cl + +debug = False + +def compile_kernels(src_path, cmd_line, *args): + global program, ctx + src = ' ' + for entry in src_path: + with open(entry, 'r') as io: + src = src + io.read() + prg = cl.Program(ctx, src).build(options=cmd_line) + program = prg + return prg + +def create_ctx(): + global ctx, queue, PYOPENCL_COMPILER_OUTPUT, devices, group_size, iscpu + PYOPENCL_COMPILER_OUTPUT=1 + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + devices = ctx.devices + group_size = devices[0].max_work_group_size + idstring = str(devices[0]) + if 'CPU' in idstring or 'cpu' in idstring: + iscpu = True + else: + iscpu = False + + return ctx, queue + +def run_kernel(names, g_range, l_range, *arguments): + global queue, program + inputbuffers = [] + for arg in arguments: + inputbuffers.append(cl.Buffer(ctx, cl.mem_flags.COPY_HOST_PTR, hostbuf = arg)) + callstrings = [] + if l_range == 0: + l_range = None + else: + l_range = tuple(l_range) + + if debug: + print("g_range: " + str(g_range) + " l_range: " +str(l_range)) + + if debug: + for arg in arguments: + print("Argumentshape:" + str(arg.shape)) + + if type(names) is list: + for name in names: + callstrings.append('program.' + name + '(queue, tuple(g_range) , l_range, *inputbuffers)') + else: + callstrings.append('program.' + names + '(queue, tuple(g_range), l_range, *inputbuffers)') + + for callstring in callstrings: + event = eval(callstring) + #event.wait() + i = 0 + for arg in arguments: + cl.enqueue_copy(queue, arg, inputbuffers[i]).wait() + inputbuffers[i].release() + i = i +1 + return 0 + +def run_kernel_ref(name, g_range, l_range, *arg_refs): + global queue, program + callstring = 'program.' + name + '(queue, None, None, *inputbuffers)' + eval(callstring) + return + +def allocate_buf(array): + return cl.Buffer(ctx, cl.mem_flags.COPY_HOST_PTR, hostbuf = array) + +def malloc(size): + return cl.Buffer(ctx, size=size) + + + + diff --git a/python/plot_2d.py b/python/plot_2d.py new file mode 100644 index 0000000..d311051 --- /dev/null +++ b/python/plot_2d.py @@ -0,0 +1,66 @@ +import matplotlib.pyplot as plt +import math + +def plot_2D(u, CSs, NODES_X, NODES_Y, NODES_Z, plot_title, field_idx1, field_idx2): + numNodes = NODES_X*NODES_Y*NODES_Z + field = u[:,field_idx1-1:field_idx2] #Matlab counts from 1 + mag = np.zeros((numNodes, 1)) + if field.shape[1] == 1: + mag[:-1] = field[1:numNodes] + else: + for i in range(numNodes): + mag[i] = np.linalg.norm(field[i, :]) + print("Warning. This is only the absolute value") + + field_cs = np.reshape(mag, (NODES_Z, NODES_Y, NODES_X)) + + for CS in CSs: + if CS[0] == 'x': + title = plot_title + ' X-Cross-Section' + v = field_cs[:,:,int(math.floor(NODES_X)/2)] + xtext = 'Y' + ytext = 'Z' + if CS[0] == 'y': + title = plot_title + ' Y-Cross-Section' + v = field_cs[:, int(math.floor(NODES_Y)/2)] + xtext = 'X' + ytext = 'Z' + if CS[0] =='z': + ytext = 'Y' + xtext = 'X' + title = plot_title + ' Z-Cross-Section' + v = field_cs[int(math.floor(NODES_Z)/2),:,:] + + a = plt.pcolor(v) + plt.colorbar(a) + plt.title(title) + plt.xlabel(xtext + " Axis (Nodes)") + plt.ylabel(ytext + " Axis (Nodes)") + plt.show() + +def plot_fieldlines_B(u, CSs, NODES_X, NODES_Y, NODES_Z, plot_title, field_idx1, field_idx2): + CS = CSs[0] + U = u[:, field_idx1-1] + V = u[:, field_idx2-1] + U = np.reshape(U, (NODES_Z, NODES_Y, NODES_X)) + V = np.reshape(V, (NODES_Z, NODES_Y, NODES_X)) + if CS[0] == 'x': + a = NODES_Y + b = NODES_Z + U = U[:,:,int(math.floor(NODES_X)/2)] + V = V[:, :, int(math.floor(NODES_X)/2)] + + if CS[0] == 'y': + a = NODES_X + b = NODES_Z + U = U[:, int(math.floor(NODES_Y)/2),:] + V = V[:, int(math.floor(NODES_Y)/2),:] + if CS[0] == 'z': + a = NODES_X + b = NODES_Y + U = U[int(math.floor(NODES_Z)/2), :, :] + V = V[int(math.floor(NODES_Z)/2), :, :] + + Y, X = np.mgrid[0:a, 0:b] + plt.streamplot(X, Y, U, V) + plt.show() diff --git a/python/prepare_vars.py b/python/prepare_vars.py new file mode 100644 index 0000000..3ce15dd --- /dev/null +++ b/python/prepare_vars.py @@ -0,0 +1,17 @@ +# from +Balance_Laws/prepare_vars.m + + +def prepare_vars(): + global I_Mesh, I_TI, I_BalanceLaws, I_Tech, I_RunOps, I_Results + I_Mesh = {'NODES_X':0, 'NODES_Y':0, 'NODES_Z':0, 'DX':0, 'DY':0, 'DZ':0, 'XMIN':0, 'XMAX':0, 'YMAX':0, 'YMIN':0, 'ZMIN':0, 'ZMAX':0} + + I_TI = {'cfl':0, 'final_time':0, 'time_integrator':'\0', 'DT':0, 'num_steps':0} + + I_Tech = {'device':0, 'REAL':'\0','optimizations':'\0','NUM_NODES_PAD':0, 'num_groups':0,'W_SIZE':0, 'g_range':0, 'l_range':0} + + I_BalanceLaws = {'hall_term':'\0', 'g_range':0, 'l_range':0, 'NUM_CONSERVED_VARS':0, 'NUM_AUXILIARY_VARS':0, 'NUM_TOTAL_VARS':0} + + I_RunOps = {'order':0, 'operator_form':'classical', 'conservation_laws':'\0', 'testcase':'\0', 'periodicx':'\0', 'periodicy':'\0', 'periodicz':'\0', 'plot_numerical_solution':'\0', 'save_fields':False, 'save_integrals_over_time':False, 'norm':'L2'} + + I_Results = {'abs_err':0, 'rel_err':0, 'field_u':0, 'runtime':0, 'kernel_runtime':0, 'error_over_time':0, 'time':0} + diff --git a/python/save_hdf5.py b/python/save_hdf5.py new file mode 100644 index 0000000..bdf0619 --- /dev/null +++ b/python/save_hdf5.py @@ -0,0 +1,40 @@ +import h5py + +def save_all_variables(array, filename, dataset_names): + NODES_X = I_Mesh['NODES_X'] + NODES_Y = I_Mesh['NODES_Y'] + NODES_Z = I_Mesh['NODES_Z'] + array = np.reshape(array, (I_Tech['NUM_NODES_PAD'], I_BalanceLaws['NUM_TOTAL_VARS'])) + f = h5py.File(filename, "w") + for i in range(len(dataset_names)): + con_var = np.reshape(array[0:NODES_X*NODES_Y*NODES_Z, i], (NODES_Z, NODES_Y, NODES_X)) + f.create_dataset(dataset_names[i], data = con_var) + + for key in I_Tech: + f.attrs[key] = I_Tech[key] + for key in I_Mesh: + f.attrs[key] = I_Mesh[key] + for key in I_TI: + f.attrs[key] = I_TI[key] + for key in I_RunOps: + f.attrs[key] = I_RunOps[key] + for key in I_BalanceLaws: + f.attrs[key] = I_BalanceLaws[key] + + f.close() + +def reload_all_variables(array, filename, dataset_names): + NODES_X = I_Mesh['NODES_X'] + NODES_Y = I_Mesh['NODES_Y'] + NODES_Z = I_Mesh['NODES_Z'] + f = h5py.File(filename, "r") + ar_shape = array.shape + ar = np.reshape(array, (I_Tech['NUM_NODES_PAD'], I_BalanceLaws['NUM_TOTAL_VARS'])) + + for i in range(len(dataset_names)): + dset = f[dataset_names[i]] + ar[0:NODES_X*NODES_Y*NODES_Z, i] = np.reshape(dset, NODES_Z * NODES_Y* NODES_Z) + + array = np.reshape(ar, ar_shape) + return array +