A commit with a series of linked changes. It was found that imposing …

…the boundary conditions that the distribution function must satisfy in the mass matrix led to numerical instability at long times. This was confirmed by removing the boundary condition imposition from the mass matrix. The addition of numerical conserving terms to the weak-form operator then provided a numerical operator that appears to have an approximate H-theorem and a stable long-time solution. Imposing numerical conserving terms alone on the old form of the mass matrix was not adequate. Here, we remove all code associated with the option impose_zero_gradient_BC, including unused functions in fokker_planck_calculus.jl. In addition to this, we also make changes in the calculus.jl and gauss_legendre.jl modules so that the comments regarding imposing boundary conditions via the 1D mass matrix are removed. We also make sure that weak-form matrices that are used only in numerical differentiation (i.e., K in M.f" = K.f, with M the mass matrix) include explicit integration-by-parts boundary condition terms so that the weak-form differentiation gives comparable or better results to the interpolation derivative method for a first and second order derivative, even when the function f has not gone to 0 by the end of the domain. This functionality can be switched by changing "K_with_BC_terms" -> "K", "L_with_BC_terms" -> "L" and toggling explicit_BC_terms = true -> explicit_BC_terms = false in the appropriate places.
mabarnes · Nov 15, 2023 · 9c8b1d3 · 9c8b1d3
1 parent 58b2494
commit 9c8b1d3
Show file tree

Hide file tree

Showing 6 changed files with 30 additions and 276 deletions.
diff --git a/2D_FEM_assembly_test.jl b/2D_FEM_assembly_test.jl
@@ -136,7 +136,7 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
 
 
     function test_weak_form_collisions(ngrid,nelement_vpa,nelement_vperp;
-        Lvpa=12.0,Lvperp=6.0,plot_test_output=false,impose_zero_gradient_BC=false,
+        Lvpa=12.0,Lvperp=6.0,plot_test_output=false,
         test_parallelism=false,test_self_operator=true,
         test_dense_construction=false,standalone=false,
         use_Maxwellian_Rosenbluth_coefficients=false,
@@ -145,7 +145,6 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
         algebraic_solve_for_d2Gdvperp2=false)
         # define inputs needed for the test
         #plot_test_output = false#true
-        #impose_zero_gradient_BC = false#true
         #test_parallelism = false#true
         #test_self_operator = true
         #test_dense_construction = false#true
@@ -198,7 +197,6 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
         KKpar2D_with_BC_terms_sparse = fkpl_arrays.KKpar2D_with_BC_terms_sparse
         KKperp2D_with_BC_terms_sparse = fkpl_arrays.KKperp2D_with_BC_terms_sparse
         lu_obj_MM = fkpl_arrays.lu_obj_MM
-        lu_obj_MMZG = fkpl_arrays.lu_obj_MMZG
         finish_init_time = now()
 
         fvpavperp = Array{mk_float,2}(undef,vpa.n,vperp.n)
@@ -233,22 +231,9 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
         # multiply by KKpar2D and fill dfc
         mul!(dfc,KKpar2D_with_BC_terms_sparse,fc)
         mul!(dgc,KKperp2D_with_BC_terms_sparse,fc)
-        if impose_zero_gradient_BC
-            # enforce zero bc  
-            enforce_zero_bc!(fc,vpa,vperp,impose_BC_at_zero_vperp=true)
-            enforce_zero_bc!(gc,vpa,vperp,impose_BC_at_zero_vperp=true)
-            # invert mass matrix and fill fc
-            fc = lu_obj_MMZG \ dfc
-            gc = lu_obj_MMZG \ dgc
-        else
-            # enforce zero bc  
-            #enforce_zero_bc!(fc,vpa,vperp,impose_BC_at_zero_vperp=false)
-            #enforce_zero_bc!(gc,vpa,vperp,impose_BC_at_zero_vperp=false)
-            # invert mass matrix and fill fc
-            fc = lu_obj_MM \ dfc
-            gc = lu_obj_MM \ dgc
-        end
-        #fc = cholesky_obj \ dfc
+        # invert mass matrix and fill fc
+        fc = lu_obj_MM \ dfc
+        gc = lu_obj_MM \ dgc
         #print_vector(fc,"fc",nc_global)
         # unravel
         ravel_c_to_vpavperp!(d2fvpavperp_dvpa2_num,fc,nc_global,vpa.n)
@@ -346,7 +331,6 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
                                              fkpl_arrays,
                                              vperp, vpa, vperp_spectral, vpa_spectral,
                                              test_assembly_serial=test_parallelism,
-                                             impose_zero_gradient_BC=impose_zero_gradient_BC,
                                              use_Maxwellian_Rosenbluth_coefficients=use_Maxwellian_Rosenbluth_coefficients,
                                              use_Maxwellian_field_particle_distribution=use_Maxwellian_field_particle_distribution,
                                              algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2)
@@ -467,7 +451,6 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
     end
 
     function run_assembly_test(; ngrid=5, nelement_list = [8],
-        impose_zero_gradient_BC= false,
         plot_scan=true,
         plot_test_output = false,
         use_Maxwellian_Rosenbluth_coefficients=false,
@@ -481,7 +464,6 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
         #ngrid = 5
         #plot_scan = true
         #plot_test_output = true#false
-        #impose_zero_gradient_BC = false
         #test_parallelism = false
         test_self_operator = true
         #test_dense_construction = false
@@ -541,7 +523,6 @@ using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary
             nelement_vperp = nelement
             fkerr, calculate_times[iscan], init_times[iscan] = test_weak_form_collisions(ngrid,nelement_vpa,nelement_vperp,
             plot_test_output=plot_test_output,
-            impose_zero_gradient_BC=impose_zero_gradient_BC,
             test_parallelism=test_parallelism,
             test_self_operator=test_self_operator,
             test_dense_construction=test_dense_construction,

diff --git a/GaussLobattoLegendre_test.jl b/GaussLobattoLegendre_test.jl
@@ -147,18 +147,14 @@ using moment_kinetics.calculus: derivative!, second_derivative!, laplacian_deriv
             println("max(d2f_err) (double first derivative by interpolation): ",maximum(d2f_err))  
             if y.name == "vpa"
                 mul!(b,y_spectral.S_matrix,f_exact)
-                #b[1] -= f_exact[1]
-                #b[y.n] += f_exact[y.n]
-                gausslegendre_mass_matrix_solve!(df_num,b,y.name,y_spectral)
+                gausslegendre_mass_matrix_solve!(df_num,b,y_spectral)
                 @. df_err = df_num - df_exact
                 #println("df_num (weak form): ",df_num)
                 #println("df_exact (weak form): ",df_exact)
                 println("max(df_err) (weak form): ",maximum(df_err))
-                #second_derivative!(d2f_num, f_exact, y, y_spectral)
-                mul!(b,y_spectral.K_matrix,f_exact)
-                #b[1] -= sum(y_spectral.lobatto.Dmat[1,:].*f_exact[1:y.ngrid])/y.element_scale[1]
-                #b[y.n] += sum(y_spectral.lobatto.Dmat[y.ngrid,:].*f_exact[y.n+1-y.ngrid:y.n])/y.element_scale[y.nelement_local]
-                gausslegendre_mass_matrix_solve!(d2f_num,b,y.name,y_spectral)
+                second_derivative!(d2f_num, f_exact, y, y_spectral)
+                #mul!(b,y_spectral.K_matrix,f_exact)
+                #gausslegendre_mass_matrix_solve!(d2f_num,b,y_spectral)
                 @. d2f_err = abs(d2f_num - d2f_exact) #(0.5*y.L/y.nelement_global)*
                 #println(d2f_num)
                 #println(d2f_exact)
@@ -170,18 +166,16 @@ using moment_kinetics.calculus: derivative!, second_derivative!, laplacian_deriv
             elseif y.name == "vperp"
                 #println("condition: ",cond(y_spectral.mass_matrix)) 
                 mul!(b,y_spectral.S_matrix,g_exact)
-                #b[y.n] += y.grid[y.n]*g_exact[y.n]    
-                gausslegendre_mass_matrix_solve!(divg_num,b,y.name,y_spectral)
+                gausslegendre_mass_matrix_solve!(divg_num,b,y_spectral)
                 @. divg_err = abs(divg_num - divg_exact)
                 #println("divg_b (weak form): ",b)
                 #println("divg_num (weak form): ",divg_num)
                 #println("divg_exact (weak form): ",divg_exact)
                 println("max(divg_err) (weak form): ",maximum(divg_err))
 
-                #second_derivative!(d2f_num, f_exact, y, y_spectral)
-                mul!(b,y_spectral.K_matrix,f_exact)
-                #b[y.n] += y.grid[y.n]*sum(y_spectral.lobatto.Dmat[y.ngrid,:].*f_exact[y.n+1-y.ngrid:y.n])/y.element_scale[y.nelement_local]
-                gausslegendre_mass_matrix_solve!(d2f_num,b,y.name,y_spectral)
+                second_derivative!(d2f_num, f_exact, y, y_spectral)
+                #mul!(b,y_spectral.K_matrix,f_exact)
+                #gausslegendre_mass_matrix_solve!(d2f_num,b,y_spectral)
                 @. d2f_err = abs(d2f_num - d2f_exact) #(0.5*y.L/y.nelement_global)*
                 #println(d2f_num)
                 #println(d2f_exact)
@@ -191,10 +185,9 @@ using moment_kinetics.calculus: derivative!, second_derivative!, laplacian_deriv
                 outfile = "vperp_second_derivative_test.pdf"
                 savefig(outfile)
 
-                mul!(b,y_spectral.L_matrix,h_exact)
-                #b[y.n] += y.grid[y.n]*sum(y_spectral.lobatto.Dmat[y.ngrid,:].*h_exact[y.n+1-y.ngrid:y.n])/y.element_scale[y.nelement_local]
-                #laplacian_derivative!(laph_num, h_exact, y, y_spectral)
-                gausslegendre_mass_matrix_solve!(laph_num,b,y.name,y_spectral)
+                laplacian_derivative!(laph_num, h_exact, y, y_spectral)
+                #mul!(b,y_spectral.L_matrix,h_exact)
+                #gausslegendre_mass_matrix_solve!(laph_num,b,y_spectral)
                 @. laph_err = abs(laph_num - laph_exact) #(0.5*y.L/y.nelement_global)*
                 #println(b[1:10])
                 #println(laph_num)

diff --git a/src/calculus.jl b/src/calculus.jl
@@ -44,7 +44,7 @@ function second_derivative!(d2f, f, coord, spectral::gausslegendre_info)
     # at element boundaries, use the average of the derivatives from neighboring elements.
     derivative_elements_to_full_grid!(coord.scratch, coord.scratch_2d, coord)
     # solve weak form problem M * d2f = K * f
-    gausslegendre_mass_matrix_solve!(d2f,coord.scratch,coord.name,spectral)
+    gausslegendre_mass_matrix_solve!(d2f,coord.scratch,spectral)
 end
 
 function laplacian_derivative!(d2f, f, coord, spectral::gausslegendre_info)
@@ -54,7 +54,7 @@ function laplacian_derivative!(d2f, f, coord, spectral::gausslegendre_info)
     # at element boundaries, use the average of the derivatives from neighboring elements.
     derivative_elements_to_full_grid!(coord.scratch, coord.scratch_2d, coord)
     # solve weak form problem M * d2f = K * f
-    gausslegendre_mass_matrix_solve!(d2f,coord.scratch,coord.name,spectral)
+    gausslegendre_mass_matrix_solve!(d2f,coord.scratch,spectral)
 end
 """
 Chebyshev transform f to get Chebyshev spectral coefficients and use them to calculate f'

diff --git a/src/fokker_planck.jl b/src/fokker_planck.jl
@@ -34,8 +34,6 @@ using ..fokker_planck_calculus: allocate_rosenbluth_potential_boundary_data
 using ..fokker_planck_calculus: fokkerplanck_arrays_struct, fokkerplanck_weakform_arrays_struct
 using ..fokker_planck_calculus: assemble_matrix_operators_dirichlet_bc
 using ..fokker_planck_calculus: assemble_matrix_operators_dirichlet_bc_sparse
-using ..fokker_planck_calculus: assemble_matrix_operators_dirichlet_bc_plus_vperp_zero_gradient
-using ..fokker_planck_calculus: assemble_matrix_operators_dirichlet_bc_plus_vperp_zero_gradient_sparse
 using ..fokker_planck_calculus: assemble_explicit_collision_operator_rhs_serial!
 using ..fokker_planck_calculus: assemble_explicit_collision_operator_rhs_parallel!
 using ..fokker_planck_calculus: assemble_explicit_collision_operator_rhs_parallel_analytical_inputs!
@@ -107,17 +105,14 @@ function init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_sp
         LP2D_sparse, LV2D_sparse, LB2D_sparse, KPperp2D_sparse,
         PUperp2D_sparse, PPparPUperp2D_sparse, PPpar2D_sparse,
         MMparMNperp2D_sparse = assemble_matrix_operators_dirichlet_bc(vpa,vperp,vpa_spectral,vperp_spectral)
-        MM2DZG_sparse = assemble_matrix_operators_dirichlet_bc_plus_vperp_zero_gradient(vpa,vperp,vpa_spectral,vperp_spectral)
     else
         MM2D_sparse, KKpar2D_sparse, KKperp2D_sparse,
         KKpar2D_with_BC_terms_sparse, KKperp2D_with_BC_terms_sparse,
         LP2D_sparse, LV2D_sparse, LB2D_sparse, KPperp2D_sparse,
         PUperp2D_sparse, PPparPUperp2D_sparse, PPpar2D_sparse,
         MMparMNperp2D_sparse = assemble_matrix_operators_dirichlet_bc_sparse(vpa,vperp,vpa_spectral,vperp_spectral)
-        MM2DZG_sparse = assemble_matrix_operators_dirichlet_bc_plus_vperp_zero_gradient_sparse(vpa,vperp,vpa_spectral,vperp_spectral)
     end
     lu_obj_MM = lu(MM2D_sparse)
-    lu_obj_MMZG = lu(MM2DZG_sparse)
     lu_obj_LP = lu(LP2D_sparse)
     lu_obj_LV = lu(LV2D_sparse)
     lu_obj_LB = lu(LB2D_sparse)
@@ -159,8 +154,8 @@ function init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_sp
     fka = fokkerplanck_weakform_arrays_struct(bwgt,rpbd,MM2D_sparse,KKpar2D_sparse,KKperp2D_sparse,
                                            KKpar2D_with_BC_terms_sparse,KKperp2D_with_BC_terms_sparse,
                                            LP2D_sparse,LV2D_sparse,LB2D_sparse,PUperp2D_sparse,PPparPUperp2D_sparse,
-                                           PPpar2D_sparse,MMparMNperp2D_sparse,KPperp2D_sparse,MM2DZG_sparse,
-                                           lu_obj_MM,lu_obj_MMZG,lu_obj_LP,lu_obj_LV,lu_obj_LB,
+                                           PPpar2D_sparse,MMparMNperp2D_sparse,KPperp2D_sparse,
+                                           lu_obj_MM,lu_obj_LP,lu_obj_LV,lu_obj_LB,
                                            YY_arrays, S_dummy, Q_dummy, rhsvpavperp, rhsc, rhqc, sc, qc,
                                            CC, HH, dHdvpa, dHdvperp, dGdvperp, d2Gdvperp2, d2Gdvpa2, d2Gdvperpdvpa,
                                            FF, dFdvpa, dFdvperp)
@@ -379,7 +374,7 @@ Function for evaluating Css' = Css'[Fs,Fs']
 function fokker_planck_collision_operator_weak_form!(ffs_in,ffsp_in,ms,msp,nussp,
                                              fkpl_arrays::fokkerplanck_weakform_arrays_struct,
                                              vperp, vpa, vperp_spectral, vpa_spectral;
-                                             test_assembly_serial=false,impose_zero_gradient_BC=false,
+                                             test_assembly_serial=false,
                                              use_Maxwellian_Rosenbluth_coefficients=false,
                                              use_Maxwellian_field_particle_distribution=false,
                                              skip_Rosenbluth_potential_calculation=false,
@@ -400,9 +395,7 @@ function fokker_planck_collision_operator_weak_form!(ffs_in,ffsp_in,ms,msp,nussp
     PPpar2D_sparse = fkpl_arrays.PPpar2D_sparse
     MMparMNperp2D_sparse = fkpl_arrays.MMparMNperp2D_sparse
     KPperp2D_sparse = fkpl_arrays.KPperp2D_sparse
-    MM2DZG_sparse = fkpl_arrays.MM2DZG_sparse
     lu_obj_MM = fkpl_arrays.lu_obj_MM
-    lu_obj_MMZG = fkpl_arrays.lu_obj_MMZG
     lu_obj_LP = fkpl_arrays.lu_obj_LP
     lu_obj_LV = fkpl_arrays.lu_obj_LV
     lu_obj_LB = fkpl_arrays.lu_obj_LB
@@ -540,15 +533,8 @@ function fokker_planck_collision_operator_weak_form!(ffs_in,ffsp_in,ms,msp,nussp
     # solve the collision operator matrix eq
     begin_serial_region()
     @serial_region begin
-        if impose_zero_gradient_BC
-            enforce_zero_bc!(rhsc,vpa,vperp,impose_BC_at_zero_vperp=true)
-            # invert mass matrix and fill fc
-            sc .= lu_obj_MMZG \ rhsc
-        else
-            #enforce_zero_bc!(rhsc,vpa,vperp)
-            # invert mass matrix and fill fc
-            sc .= lu_obj_MM \ rhsc
-        end
+        # invert mass matrix and fill fc
+        sc .= lu_obj_MM \ rhsc
     end
     ravel_c_to_vpavperp_parallel!(CC,sc,vpa.n)
     return nothing