Skip to content

Commit

Permalink
continue wall bc Jacobian
Browse files Browse the repository at this point in the history 
start adding effect of delta_ppar. some other small fixes
  • Loading branch information
@johnomotani
johnomotani committed Nov 25, 2024
1 parent 59a469a commit 8e826b1
Show file tree
Hide file tree
Showing 2 changed files with 165 additions and 22 deletions.
148 changes: 133 additions & 15 deletions moment_kinetics/src/electron_kinetic_equation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2942,18 +2942,21 @@ end
end

"""
add_wall_boundary_condition_to_Jacobian!(jacobian, phi, vthe, upar, z, vperp, vpa,
vperp_spectral, vpa_spectral, vpa_adv,
moments, vpa_diffusion, me_over_mi, ir)
add_wall_boundary_condition_to_Jacobian!(jacobian, phi, pdf, ppar, vthe, upar, z,
vperp, vpa, vperp_spectral, vpa_spectral,
vpa_adv, moments, vpa_diffusion, me_over_mi,
ir)
"""
@timeit global_timer add_wall_boundary_condition_to_Jacobian!(
jacobian, phi, vthe, upar, z, vperp, vpa, vperp_spectral,
vpa_spectral, vpa_adv, moments, vpa_diffusion, me_over_mi,
ir) = begin
jacobian, phi, pdf, ppar, vthe, upar, z, vperp, vpa,
vperp_spectral, vpa_spectral, vpa_adv, moments, vpa_diffusion,
me_over_mi, ir) = begin
if z.bc != "wall"
return nothing
end

pdf_size = z.n * vperp.n * vpa.n

if z.irank == 0
begin_vperp_region()
@loop_vperp ivperp begin
Expand All @@ -2968,6 +2971,7 @@ end

jac_range = (ivperp-1)*vpa.n+1 : ivperp*vpa.n
jacobian_zbegin = @view jacobian[jac_range,jac_range]
jacobian_zbegin_ppar = @view jacobian[jac_range,pdf_size+1]

vpa_unnorm, u_over_vt, vcut, minus_vcut_ind, sigma, sigma_ind, sigma_fraction,
element_with_zero, element_with_zero_boundary, last_point_near_zero,
Expand All @@ -2977,7 +2981,7 @@ end
plus_vcut_ind = searchsortedlast(vpa_unnorm, vcut)
# vcut_fraction is the fraction of the distance between plus_vcut_ind and
# plus_vcut_ind+1 where vcut is.
vcut_fraction = (vcut - vpa_unnorm[plus_vcut_ind]) / (vpa_unnorm[plus_vcut_ind+1] - vpa_unnorm[plus_vcut_ind])
vcut_fraction = get_plus_vcut_fraction(vcut, plus_vcut_ind, vpa_unnorm)
if vcut_fraction > 0.5
last_nonzero_ind = plus_vcut_ind + 1
else
Expand All @@ -3001,12 +3005,68 @@ end
else
jacobian_zbegin[last_nonzero_ind,:] .*= vcut_fraction + 0.5
end

# Fill in elements giving response to changes in electron_ppar
# A change to ppar results in a shift in the location of w_∥(v_∥=0).
# The interpolated values of g_e(w_∥) that are filled by the boundary
# condition are (dropping _∥ subscripts for the remaining comments in this
# if-clause)
# g(w|v>0) = g(2*u/vth - w)
# So
# δg(w|v>0) = dg/dw(2*u/vth - w) * (-2*u/vth^2) * δvth
# = dg/dw(2*u/vth - w) * (-2*u/vth^2) * δ(sqrt(2*ppar/n)
# = dg/dw(2*u/vth - w) * (-2*u/vth^2) * δppar * vth / (2*ppar)
# = -u/vth/ppar * dg/dw(2*u/vth - w) * δppar

dpdfdw_near_zero = @view vpa.scratch[sigma_ind:last_point_near_zero]
@views interpolate_symmetric!(dpdfdw_near_zero,
vpa_unnorm[sigma_ind:last_point_near_zero],
pdf[element_with_zero_boundary:sigma_ind-1,1,1,ir],
vpa_unnorm[element_with_zero_boundary:sigma_ind-1],
Val(1))
@. jacobian_zbegin_ppar[sigma_ind:last_point_near_zero] =
-upar[1] / vthe[1] / ppar[1] * dpdfdw_near_zero

reversed_dpdfdw_far_from_zero = @view vpa.scratch[last_point_near_zero+1:last_nonzero_ind]
@views interpolate_to_grid_1d!(reversed_dpdfdw_far_from_zero,
reversed_wpa_of_minus_vpa[vpa.n-last_nonzero_ind+1:vpa.n-last_point_near_zero],
pdf[:,1,1,ir], vpa, vpa_spectral, Val(1))
reverse!(reversed_dpdfdw_far_from_zero)
# Note that because we calculated the derivative of the interpolating
# function, and then reversed the results, we need to multiply the derivative
# by -1.
@. jacobian_zbegin_ppar[last_point_near_zero+1:last_nonzero_ind] =
upar[1] / vthe[1] / ppar[1] * reversed_dpdfdw_far_from_zero

# The change in electron_ppar also changes the position of
# wcut = (vcut - upar)/vthe
# as vcut does not change (within the Krylov iteration where this
# preconditioner matrix is used), but vthe does because of the change in
# electron_ppar. We actually use vcut_fraction calculated from vpa_unnorm, so
# it is most convenient to consider:
# v = vthe * w + upar
# δv = δvthe * w
# = δvthe * (v - upar)/vthe
# = δppar * vthe / (2*ppar) * (v - upar)/vthe
# = δppar * (v - upar) / 2 / ppar
# with vl and vu the values of v at the grid points below and above vcut
# vcut_fraction = (vcut - vl) / (vu - vl)
# δvcut_fraction = -(vcut - vl) / (vu - vl)^2 * (δvu - δvl) - δvl / (vu - vl)
# δvcut_fraction = [-(vcut - vl) / (vu - vl)^2 * (vu - vl) - (vl - upar) / (vu - vl)] * δppar / 2 / ppar
# δvcut_fraction = [-(vcut - vl) / (vu - vl) - (vl - upar) / (vu - vl)] * δppar / 2 / ppar
# δvcut_fraction = -(vcut - upar) / (vu - vl) / 2 / ppar * δppar
interpolated_pdf_at_last_nonzero_ind = @view vpa.scratch[1:1]
@views interpolate_to_grid_1d!(interpolated_pdf_at_last_nonzero_ind,
vpa_unnorm[last_nonzero_ind:last_nonzero_ind],
pdf[:,1,1,ir], vpa, vpa_spectral)

delta_vcut_fraction_over_delta_ppar = -(vcut - upar[1]) / (vpa_unnorm[plus_vcut_ind+1] - vpa_unnorm[plus_vcut_ind]) / 2.0 / ppar[1]
jacobian_zbegin_ppar[last_nonzero_ind] += interpolated_pdf_at_last_nonzero_ind[] * delta_vcut_fraction_over_delta_ppar
end
end

if z.irank == z.nrank - 1
begin_vperp_region()
pdf_size = z.n * vperp.n * vpa.n
@loop_vperp ivperp begin
# Skip vperp boundary points.
if vperp.n > 1 && ivperp == vperp.n
Expand All @@ -3019,18 +3079,19 @@ end

jac_range = pdf_size-vperp.n*vpa.n+(ivperp-1)*vpa.n+1 : pdf_size-vperp.n*vpa.n+ivperp*vpa.n
jacobian_zend = @view jacobian[jac_range,jac_range]
jacobian_zend_ppar = @view jacobian[jac_range,pdf_size+z.n]

vpa_unnorm, u_over_vt, vcut, plus_vcut_ind, sigma, sigma_ind, sigma_fraction,
element_with_zero, element_with_zero_boundary, first_point_near_zero,
reversed_wpa_of_minus_vpa = get_cutoff_params_upper(upar, vthe, phi,
me_over_mi, vpa, ir)

minus_vcut_ind = searchsortedfirst(vpa_unnorm, -vcut)
# vcut_fraction is the fraction of the distance between minus_vcut_ind and
# minus_vcut_ind-1 where -vcut is.
vcut_fraction = (-vcut - vpa_unnorm[minus_vcut_ind]) / (vpa_unnorm[minus_vcut_ind-1] - vpa_unnorm[minus_vcut_ind])
# vcut_fraction is the fraction of the distance between minus_vcut_ind-1 and
# minus_vcut_ind where -vcut is.
vcut_fraction = get_minus_vcut_fraction(vcut, minus_vcut_ind, vpa_unnorm)

if vcut_fraction > 0.5
if vcut_fraction < 0.5
first_nonzero_ind = minus_vcut_ind - 1
else
first_nonzero_ind = minus_vcut_ind
Expand All @@ -3048,11 +3109,68 @@ end
reversed_wpa_of_minus_vpa[vpa.n-first_point_near_zero+2:vpa.n-first_nonzero_ind+1],
vpa, vpa_spectral)

if vcut_fraction > 0.5
jacobian_zend[first_nonzero_ind,:] .*= vcut_fraction - 0.5
if vcut_fraction < 0.5
jacobian_zend[first_nonzero_ind,:] .*= 0.5 - vcut_fraction
else
jacobian_zend[first_nonzero_ind,:] .*= vcut_fraction + 0.5
jacobian_zend[first_nonzero_ind,:] .*= 1.5 - vcut_fraction
end

# Fill in elements giving response to changes in electron_ppar
# A change to ppar results in a shift in the location of w_∥(v_∥=0).
# The interpolated values of g_e(w_∥) that are filled by the boundary
# condition are (dropping _∥ subscripts for the remaining comments in this
# if-clause)
# g(w|v<0) = g(2*u/vth - w)
# So
# δg(w|v<0) = dg/dw(2*u/vth - w) * (-2*u/vth^2) * δvth
# = dg/dw(2*u/vth - w) * (-2*u/vth^2) * δ(sqrt(2*ppar/n)
# = dg/dw(2*u/vth - w) * (-2*u/vth^2) * δppar * vth / (2*ppar)
# = -u/vth/ppar * dg/dw(2*u/vth - w) * δppar

dpdfdw_near_zero = @view vpa.scratch[first_point_near_zero:sigma_ind]
@views interpolate_symmetric!(dpdfdw_near_zero,
vpa_unnorm[first_point_near_zero:sigma_ind],
pdf[sigma_ind+1:element_with_zero_boundary,1,end,ir],
vpa_unnorm[sigma_ind+1:element_with_zero_boundary],
Val(1))
@. jacobian_zend_ppar[first_point_near_zero:sigma_ind] =
-upar[end] / vthe[end] / ppar[end] * dpdfdw_near_zero

reversed_dpdfdw_far_from_zero = @view vpa.scratch[first_nonzero_ind:first_point_near_zero-1]
@views interpolate_to_grid_1d!(reversed_dpdfdw_far_from_zero,
reversed_wpa_of_minus_vpa[vpa.n-first_point_near_zero+2:vpa.n-first_nonzero_ind+1],
pdf[:,1,end,ir], vpa, vpa_spectral, Val(1))
reverse!(reversed_dpdfdw_far_from_zero)
# Note that because we calculated the derivative of the interpolating
# function, and then reversed the results, we need to multiply the derivative
# by -1.
@. jacobian_zend_ppar[first_nonzero_ind:first_point_near_zero-1] =
upar[end] / vthe[end] / ppar[end] * reversed_dpdfdw_far_from_zero

# The change in electron_ppar also changes the position of
# wcut = (vcut - upar)/vthe
# as vcut does not change (within the Krylov iteration where this
# preconditioner matrix is used), but vthe does because of the change in
# electron_ppar. We actually use vcut_fraction calculated from vpa_unnorm, so
# it is most convenient to consider:
# v = vthe * w + upar
# δv = δvthe * w
# = δvthe * (v - upar)/vthe
# = δppar * vthe / (2*ppar) * (v - upar)/vthe
# = δppar * (v - upar) / 2 / ppar
# with vl and vu the values of v at the grid points below and above vcut
# vcut_fraction = (vcut - vl) / (vu - vl)
# δvcut_fraction = -(vcut - vl) / (vu - vl)^2 * (δvu - δvl) - δvl / (vu - vl)
# δvcut_fraction = [-(vcut - vl) / (vu - vl)^2 * (vu - vl) - (vl - upar) / (vu - vl)] * δppar / 2 / ppar
# δvcut_fraction = [-(vcut - vl) / (vu - vl) - (vl - upar) / (vu - vl)] * δppar / 2 / ppar
# δvcut_fraction = -(vcut - upar) / (vu - vl) / 2 / ppar * δppar
interpolated_pdf_at_first_nonzero_ind = @view vpa.scratch[1:1]
@views interpolate_to_grid_1d!(interpolated_pdf_at_first_nonzero_ind,
vpa_unnorm[first_nonzero_ind:first_nonzero_ind],
pdf[:,1,1,ir], vpa, vpa_spectral)

delta_vcut_fraction_over_delta_ppar = -(vcut - upar[end]) / (vpa_unnorm[minus_vcut_ind] - vpa_unnorm[minus_vcut_ind-1]) / 2.0 / ppar[1]
jacobian_zend_ppar[first_nonzero_ind] += interpolated_pdf_at_first_nonzero_ind[] * delta_vcut_fraction_over_delta_ppar
end
end
end
Expand Down
39 changes: 32 additions & 7 deletions moment_kinetics/test/jacobian_matrix_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ using moment_kinetics.moment_constraints: electron_implicit_constraint_forcing!,
add_electron_implicit_constraint_forcing_to_z_only_Jacobian!,
add_electron_implicit_constraint_forcing_to_v_only_Jacobian!,
hard_force_moment_constraints!
using moment_kinetics.timer_utils: reset_mk_timers!
using moment_kinetics.type_definitions: mk_float
using moment_kinetics.velocity_moments: calculate_electron_moment_derivatives_no_r!,
integrate_over_vspace
Expand Down Expand Up @@ -210,6 +211,9 @@ test_input = OptionsDict("output" => OptionsDict("run_name" => "jacobian_matrix"
)

function get_mk_state(test_input)
# Reset timers in case there was a previous run which did not clean them up.
reset_mk_timers!()

mk_state = nothing
quietoutput() do
mk_state = setup_moment_kinetics(test_input; skip_electron_solve=true)
Expand Down Expand Up @@ -3620,8 +3624,6 @@ function test_electron_wall_bc(test_input; atol=(5.0*epsilon)^2)
f = @view pdf.electron.norm[:,:,:,ir]
begin_serial_region()
@serial_region begin
@. delta_p = p_amplitude * sin(2.0*π*test_wavenumber*z.grid/z.L)

# Make sure initial condition has some z-variation. As f is 'moment kinetic' this
# means f must have a non-Maxwellian part that varies in z.
f .*= 1.0 .+ 1.0e-4 .* reshape(vpa.grid.^3, vpa.n, 1, 1) .* reshape(sin.(2.0.*π.*z.grid./z.L), 1, 1, z.n)
Expand All @@ -3636,6 +3638,8 @@ function test_electron_wall_bc(test_input; atol=(5.0*epsilon)^2)
# at the z-boundaries
begin_serial_region()
@serial_region begin
@. delta_p = p_amplitude

delta_f .= f_amplitude .*
reshape(exp.(sin.(2.0.*π.*test_wavenumber.*vpa.grid./vpa.L)) .- 1.0, vpa.n, 1, 1) .*
f
Expand Down Expand Up @@ -3678,9 +3682,9 @@ function test_electron_wall_bc(test_input; atol=(5.0*epsilon)^2)
end

add_wall_boundary_condition_to_Jacobian!(
jacobian_matrix, phi, vth, upar, z, vperp, vpa, vperp_spectral, vpa_spectral,
vpa_advect, moments, num_diss_params.electron.vpa_dissipation_coefficient, me,
ir)
jacobian_matrix, phi, f, ppar, vth, upar, z, vperp, vpa, vperp_spectral,
vpa_spectral, vpa_advect, moments,
num_diss_params.electron.vpa_dissipation_coefficient, me, ir)

# # Test 'ADI Jacobians' before other tests, because residual_func() may modify some
# # variables (vth, etc.).
Expand Down Expand Up @@ -3764,8 +3768,7 @@ function test_electron_wall_bc(test_input; atol=(5.0*epsilon)^2)
before = copy(this_f)
# enforce the boundary condition(s) on the electron pdf
@views enforce_boundary_condition_on_electron_pdf!(
this_f, phi, moments.electron.vth[:,ir],
moments.electron.upar[:,ir], z, vperp, vpa, vperp_spectral,
this_f, phi, vth, upar, z, vperp, vpa, vperp_spectral,
vpa_spectral, vpa_advect, moments,
num_diss_params.electron.vpa_dissipation_coefficient > 0.0,
composition.me_over_mi, ir; bc_constraints=false,
Expand Down Expand Up @@ -3837,7 +3840,29 @@ before = copy(this_f)
end
end

@testset "δp only" begin
residual_func!(original_residual, f, ppar)
residual_func!(perturbed_residual, copy(f), ppar.+delta_p)

begin_serial_region()
@serial_region begin
delta_state = zeros(mk_float, total_size)
delta_state[pdf_size+1:end] .= vec(delta_p)
residual_update_with_Jacobian = jacobian_matrix * delta_state
perturbed_with_Jacobian = vec(original_residual) .+ residual_update_with_Jacobian[1:pdf_size]

# Check ppar did not get perturbed by the Jacobian
@test elementwise_isapprox(residual_update_with_Jacobian[pdf_size+1:end],
zeros(p_size); atol=1.0e-15)

# Use an absolute tolerance for this test because if we used a norm_factor
# like the other tests, it would be zero to machine precision at some
# points.
@test elementwise_isapprox(perturbed_residual,
reshape(perturbed_with_Jacobian, vpa.n, vperp.n, z.n);
rtol=0.0, atol=atol)
end
end
# @testset "δp only" begin
# residual_func!(original_residual, f, ppar)
# residual_func!(perturbed_residual, f, ppar .+ delta_p)
Expand Down

0 comments on commit 8e826b1

Please sign in to comment.