Merge pull request #224 from mabarnes/fix_collision_operator_scripts

Improve collision operator scripts to save HDF5 format error data for easier comparisons

examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml

@@ -49,6 +49,7 @@ vperp_ngrid = 5
 vperp_nelement = 4
 vperp_L = 3.0
 vperp_discretization = "gausslegendre_pseudospectral"
+vperp_bc = "zero"
 # Fokker-Planck operator requires the "gausslegendre_pseudospectral
 # options for the vpa and vperp grids
 

moment_kinetics/src/fokker_planck_test.jl

@@ -15,10 +15,11 @@ export H_Maxwellian, G_Maxwellian
 export Cssp_fully_expanded_form, calculate_collisional_fluxes
 
 export print_test_data, fkpl_error_data, allocate_error_data
-
+export save_fkpl_error_data, save_fkpl_integration_error_data
 #using Plots
 #using LaTeXStrings
-using Measures
+#using Measures
+using HDF5
 using ..type_definitions: mk_float, mk_int
 using SpecialFunctions: erf
 using ..velocity_moments: get_density
@@ -339,4 +340,76 @@ function allocate_error_data()
         moments)
 end
 
+function save_fkpl_error_data(outdir,ncore,ngrid,nelement_list,
+    max_C_err, max_H_err, max_G_err, max_dHdvpa_err, max_dHdvperp_err,
+    max_d2Gdvperp2_err, max_d2Gdvpa2_err, max_d2Gdvperpdvpa_err, max_dGdvperp_err, 
+    L2_C_err, L2_H_err, L2_G_err, L2_dHdvpa_err, L2_dHdvperp_err, L2_d2Gdvperp2_err,
+    L2_d2Gdvpa2_err, L2_d2Gdvperpdvpa_err, L2_dGdvperp_err,
+    n_err, u_err, p_err, calculate_times, init_times, expected_t_2, expected_t_3,
+    expected_diff, expected_integral)
+    filename = outdir*"fkpl_error_data_ngrid_"*string(ngrid)*"_ncore_"*string(ncore)*".h5"
+    fid = h5open(filename, "w")
+    fid["ncore"] = ncore
+    fid["ngrid"] = ngrid
+    fid["nelement_list"] = nelement_list
+    fid["max_C_err"] = max_C_err
+    fid["max_H_err"] = max_H_err
+    fid["max_G_err"] = max_G_err
+    fid["max_dHdvpa_err"] = max_dHdvpa_err
+    fid["max_dHdvperp_err"] = max_dHdvperp_err
+    fid["max_d2Gdvperp2_err"] = max_d2Gdvperp2_err
+    fid["max_d2Gdvpa2_err"] = max_d2Gdvpa2_err
+    fid["max_d2Gdvperpdvpa_err"] = max_d2Gdvperpdvpa_err
+    fid["max_dGdvperp_err"] = max_dGdvperp_err
+    fid["L2_C_err"] = L2_C_err
+    fid["L2_H_err"] = L2_H_err
+    fid["L2_G_err"] = L2_G_err
+    fid["L2_dHdvpa_err"] = L2_dHdvpa_err
+    fid["L2_dHdvperp_err"] = L2_dHdvperp_err
+    fid["L2_d2Gdvperp2_err"] = L2_d2Gdvperp2_err
+    fid["L2_d2Gdvpa2_err"] = L2_d2Gdvpa2_err
+    fid["L2_d2Gdvperpdvpa_err"] = L2_d2Gdvperpdvpa_err
+    fid["L2_dGdvperp_err"] = L2_dGdvperp_err
+    fid["n_err"] = n_err
+    fid["u_err"] = u_err
+    fid["p_err"] = p_err
+    fid["calculate_times"] = calculate_times
+    fid["init_times"] = init_times
+    fid["expected_t_2"] = expected_t_2
+    fid["expected_t_3"] = expected_t_3
+    fid["expected_diff"] = expected_diff
+    fid["expected_integral"] = expected_integral
+    close(fid)
+    println("Saving error data: ",filename)
+    return nothing
+end
+
+function save_fkpl_integration_error_data(outdir,ncore,ngrid,nelement_list,
+    max_dfsdvpa_err, max_dfsdvperp_err, max_d2fsdvperpdvpa_err,
+    max_H_err, max_G_err, max_dHdvpa_err, max_dHdvperp_err,
+    max_d2Gdvperp2_err, max_d2Gdvpa2_err, max_d2Gdvperpdvpa_err, max_dGdvperp_err, 
+    expected_diff, expected_integral)
+    filename = outdir*"fkpl_integration_error_data_ngrid_"*string(ngrid)*"_ncore_"*string(ncore)*".h5"
+    fid = h5open(filename, "w")
+    fid["ncore"] = ncore
+    fid["ngrid"] = ngrid
+    fid["nelement_list"] = nelement_list
+    fid["max_dfsdvpa_err"] = max_dfsdvpa_err
+    fid["max_dfsdvperp_err"] = max_dfsdvperp_err
+    fid["max_d2fsdvperpdvpa_err"] = max_d2fsdvperpdvpa_err
+    fid["max_H_err"] = max_H_err
+    fid["max_G_err"] = max_G_err
+    fid["max_dHdvpa_err"] = max_dHdvpa_err
+    fid["max_dHdvperp_err"] = max_dHdvperp_err
+    fid["max_d2Gdvperp2_err"] = max_d2Gdvperp2_err
+    fid["max_d2Gdvpa2_err"] = max_d2Gdvpa2_err
+    fid["max_d2Gdvperpdvpa_err"] = max_d2Gdvperpdvpa_err
+    fid["max_dGdvperp_err"] = max_dGdvperp_err
+    fid["expected_diff"] = expected_diff
+    fid["expected_integral"] = expected_integral
+    close(fid)
+    println("Saving error data: ",filename)
+    return nothing
+end
+
 end

test_scripts/2D_FEM_assembly_test.jl

@@ -28,6 +28,7 @@ using moment_kinetics.fokker_planck_test: d2Gdvpa2_Maxwellian, d2Gdvperp2_Maxwel
 using moment_kinetics.fokker_planck_test: dHdvperp_Maxwellian, dHdvpa_Maxwellian
 using moment_kinetics.fokker_planck_test: Cssp_Maxwellian_inputs
 using moment_kinetics.fokker_planck_test: print_test_data, fkpl_error_data, allocate_error_data
+using moment_kinetics.fokker_planck_test: save_fkpl_error_data
 
 using moment_kinetics.fokker_planck_calculus: elliptic_solve!
 using moment_kinetics.fokker_planck_calculus: enforce_zero_bc!, allocate_rosenbluth_potential_boundary_data
@@ -382,6 +383,7 @@ end
 
     function run_assembly_test(; ngrid=5, nelement_list = [8],
         plot_scan=true,
+        save_HDF5 = true,
         plot_test_output = false,
         use_Maxwellian_Rosenbluth_coefficients=false,
         use_Maxwellian_field_particle_distribution=false,
@@ -581,6 +583,18 @@ end
             savefig(outfile)
             println(outfile)
             println([calculate_times, init_times, expected_t_2, expected_t_3])
+
+        end
+        if global_rank[]==0 && save_HDF5
+            outdir = ""
+            ncore = global_size[]
+            save_fkpl_error_data(outdir,ncore,ngrid,nelement_list,
+                max_C_err, max_H_err, max_G_err, max_dHdvpa_err, max_dHdvperp_err,
+                max_d2Gdvperp2_err, max_d2Gdvpa2_err, max_d2Gdvperpdvpa_err, max_dGdvperp_err, 
+                L2_C_err, L2_H_err, L2_G_err, L2_dHdvpa_err, L2_dHdvperp_err, L2_d2Gdvperp2_err,
+                L2_d2Gdvpa2_err, L2_d2Gdvperpdvpa_err, L2_dGdvperp_err,
+                n_err, u_err, p_err, calculate_times, init_times, expected_t_2, expected_t_3,
+                expected, expected_integral)
         end
         finalize_comms!()
     return nothing

test_scripts/fkpl_direct_integration_test.jl

@@ -16,6 +16,7 @@ using moment_kinetics.fokker_planck_test: d2Gdvpa2_Maxwellian, dGdvperp_Maxwelli
 using moment_kinetics.fokker_planck_test: dHdvpa_Maxwellian, dHdvperp_Maxwellian, H_Maxwellian, G_Maxwellian
 using moment_kinetics.fokker_planck_test: F_Maxwellian, dFdvpa_Maxwellian, dFdvperp_Maxwellian
 using moment_kinetics.fokker_planck_test: d2Fdvpa2_Maxwellian, d2Fdvperpdvpa_Maxwellian, d2Fdvperp2_Maxwellian
+using moment_kinetics.fokker_planck_test: save_fkpl_integration_error_data
 using moment_kinetics.type_definitions: mk_float, mk_int
 using moment_kinetics.calculus: derivative!
 using moment_kinetics.velocity_moments: integrate_over_vspace, get_pressure
@@ -99,16 +100,16 @@ test_Lagrange_integral = false #true
 test_Lagrange_integral_scan = true
 
 function test_Lagrange_Rosenbluth_potentials(ngrid,nelement; standalone=true)
-    # set up grids for input Maxwellian
-    vpa, vperp, vpa_spectral, vperp_spectral =  init_grids(nelement,ngrid)
-    # set up necessary inputs for collision operator functions 
-    nvperp = vperp.n
-    nvpa = vpa.n
     # Set up MPI
     if standalone
         initialize_comms!()
     end
     setup_distributed_memory_MPI(1,1,1,1)
+    # set up grids for input Maxwellian
+    vpa, vperp, vpa_spectral, vperp_spectral =  init_grids(nelement,ngrid)
+    # set up necessary inputs for collision operator functions 
+    nvperp = vperp.n
+    nvpa = vpa.n
     looping.setup_loop_ranges!(block_rank[], block_size[];
                                    s=1, sn=1,
                                    r=1, z=1, vperp=vperp.n, vpa=vpa.n,
@@ -344,7 +345,7 @@ function test_Lagrange_Rosenbluth_potentials(ngrid,nelement; standalone=true)
     return results 
 end
 
-function test_rosenbluth_potentials_direct_integration(;ngrid=5,nelement_list=[2],plot_scan=true)
+function test_rosenbluth_potentials_direct_integration(;ngrid=5,nelement_list=[2],plot_scan=true,save_HDF5=true)
     if size(nelement_list,1) == 1
         nelement = nelement_list[1]
         test_Lagrange_Rosenbluth_potentials(ngrid,nelement,standalone=true)
@@ -427,6 +428,15 @@ function test_rosenbluth_potentials_direct_integration(;ngrid=5,nelement_list=[2
             savefig(outfile)
             println(outfile)
         end
+        if global_rank[]==0 && save_HDF5
+            outdir = ""
+            ncore = global_size[]
+            save_fkpl_integration_error_data(outdir, ncore, ngrid, nelement_list,
+                max_dfsdvpa_err, max_dfsdvperp_err, max_d2fsdvperpdvpa_err,
+                max_H_err, max_G_err, max_dHdvpa_err, max_dHdvperp_err,
+                max_d2Gdvperp2_err, max_d2Gdvpa2_err, max_d2Gdvperpdvpa_err, max_dGdvperp_err, 
+                expected, expected_integral)
+        end
         finalize_comms!()
     end
 end