Skip to content

Commit

Permalink
Changes suggested by @johnomotani to remove boundary conditions from …
Browse files Browse the repository at this point in the history 
…mass matrix inversion when carrying out differentiation, and including explicit boundary terms. Here this is handled with explicit boundary terms in the test script GaussLobattoLegendre_test.jl for 1D differentiation operations only using assembled matrices (derivative functions not yet supported). Subject to testing, these changes will be implemented across the gauss_legendre.jl and fokker_planck*.jl modules.
  • Loading branch information
@mrhardman
mrhardman committed Nov 11, 2023
1 parent f09aca8 commit 1360367
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 14 deletions.
26 changes: 16 additions & 10 deletions GaussLobattoLegendre_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,14 +147,18 @@ 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)
@. 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)
#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)
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)
@. d2f_err = abs(d2f_num - d2f_exact) #(0.5*y.L/y.nelement_global)*
#println(d2f_num)
#println(d2f_exact)
Expand All @@ -165,18 +169,19 @@ using moment_kinetics.calculus: derivative!, second_derivative!, laplacian_deriv

elseif y.name == "vperp"
#println("condition: ",cond(y_spectral.mass_matrix))
b = Array{Float64,1}(undef,y.n)
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)
@. 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)
#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)
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)
@. d2f_err = abs(d2f_num - d2f_exact) #(0.5*y.L/y.nelement_global)*
#println(d2f_num)
#println(d2f_exact)
Expand All @@ -186,9 +191,10 @@ 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)
laplacian_derivative!(laph_num, h_exact, y, y_spectral)
#gausslegendre_mass_matrix_solve!(laph_num,b,y_spectral)
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)
@. laph_err = abs(laph_num - laph_exact) #(0.5*y.L/y.nelement_global)*
#println(b[1:10])
#println(laph_num)
Expand Down
8 changes: 4 additions & 4 deletions src/gauss_legendre.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -343,11 +343,11 @@ function gausslegendre_mass_matrix_solve!(f,b,coord_name,spectral)
if coord_name == "vperp"
# enforce zero (value or gradient) boundary conditions
#b[1] = 0.0 # uncomment if bc is imposed at vperp = 0 in mass matrix
b[end] = 0.0
#b[end] = 0.0
else
# enforce zero (value or gradient) boundary conditions
b[1] = 0.0
b[end] = 0.0
#b[1] = 0.0
#b[end] = 0.0
end
# invert mass matrix system
y = spectral.mass_matrix_lu \ b
Expand Down Expand Up @@ -846,7 +846,7 @@ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2},
imin = coord.imin
imax = coord.imax
@. QQ_global = 0.0
mass_matrix = option == "M"
mass_matrix = (option == "M") && false
if coord.name == "vperp"
zero_bc_upper_boundary = true && mass_matrix
zero_bc_lower_boundary = false && mass_matrix
Expand Down

2 comments on commit 1360367

@mrhardman
Copy link
Collaborator Author

@mrhardman mrhardman commented on 1360367 Nov 11, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@johnomotani At your convenience, could you please check to see if there are still problems with the 1D differentiation test in GaussLobattoLegendre_test.jl in the intermediate integration by parts step? A brief check now suggests that the errors are much smaller (and comparable between weak and interpolation methods) for L_in = O(1). If you confirm that things are now looking OK then I will roll the changes out systematically to the other affected modules.

UPDATE: Looking at small ngrid (2,3,4) and small L_in (<1), the errors on the weak-form second derivatives still seem potentially large even for large values of nelement.

@johnomotani
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mrhardman nice! That looks like errors are now small at the boundary even when the test function has non-zero derivative at the boundary. Interestingly, there is a slightly bigger error at the boundary point than when testing with zero-ed out boundary terms and a test function that has exactly zero derivative at the boundary - I guess the numerical derivative that effectively feeds into the boundary terms introduces some error?? Anyway, it's a small error now, so fingers crossed it will be OK. Would be nice if this would make the long-time stability of your simulations better!

Please sign in to comment.