Commit 677e49c5 authored by Jannis Teunissen's avatar Jannis Teunissen

Fix sign in gradient of lsf solution

parent 0772b5a6
......@@ -33,6 +33,7 @@ module m_af_multigrid
public :: mg_box_gsrb_lpl
public :: mg_box_corr_lpl
public :: mg_box_rstr_lpl
public :: mg_box_lpl_gradient
public :: mg_sides_rb
! Methods for Laplacian with jump in coefficient between boxes
......@@ -41,6 +42,12 @@ module m_af_multigrid
public :: mg_box_corr_lpld
public :: mg_box_lpld_stencil
! Methods for Laplacian with level set function
public :: mg_box_lpllsf
public :: mg_box_lpllsf_stencil
public :: mg_box_gsrb_lpllsf
public :: mg_box_lpllsf_gradient
! To adjust operator stencils near boundaries
public :: mg_stencil_handle_boundaries
......@@ -2052,11 +2059,12 @@ contains
associate (box => tree%boxes(id))
select case(box%tag)
case (mg_normal_box, mg_ceps_box)
call mg_box_gradient(box, mg, i_fc, fac)
call mg_box_lpl_gradient(box, mg, i_fc, fac)
case (mg_lsf_box)
call mg_box_lpllsf_gradient(box, mg, i_fc, fac)
case (mg_veps_box)
error stop "Use dielectric_correct_field_fc instead"
! Should call dielectric_correct_field_fc afterwards
call mg_box_lpl_gradient(box, mg, i_fc, fac)
case (af_init_tag)
error stop "mg_auto_op: box tag not set"
case default
......@@ -2074,7 +2082,7 @@ contains
end subroutine mg_compute_phi_gradient
!> Compute the gradient of the potential
subroutine mg_box_gradient(box, mg, i_fc, fac)
subroutine mg_box_lpl_gradient(box, mg, i_fc, fac)
type(box_t), intent(inout) :: box
type(mg_t), intent(in) :: mg
integer, intent(in) :: i_fc !< Face-centered indices
......@@ -2105,7 +2113,7 @@ contains
(box%cc(1:nc, 1:nc, 1:nc+1, i_phi) - &
box%cc(1:nc, 1:nc, 0:nc, i_phi))
#endif
end subroutine mg_box_gradient
end subroutine mg_box_lpl_gradient
subroutine mg_box_gradient_norm(box, i_fc, i_norm)
type(box_t), intent(inout) :: box
......@@ -2167,7 +2175,7 @@ contains
call lsf_dist_val(v_a(1), v_b(1), v_b(2), &
mg%lsf_boundary_value, dd, val)
box%fc(IJK, 1, i_fc) = grad_sign * fac * (val - v_a(2)) / &
box%fc(IJK, 1, i_fc) = grad_sign * fac * (v_a(2) - val) / &
(box%dr(1) * dd)
end do
#elif NDIM == 2
......@@ -2185,7 +2193,7 @@ contains
call lsf_dist_val(v_a(1), v_b(1), v_b(2), &
mg%lsf_boundary_value, dd, val)
box%fc(IJK, 1, i_fc) = grad_sign * fac * (val - v_a(2)) / &
box%fc(IJK, 1, i_fc) = grad_sign * fac * (v_a(2) - val) / &
(box%dr(1) * dd)
if (box%cc(j, i, i_lsf) > 0) then
......@@ -2200,7 +2208,7 @@ contains
call lsf_dist_val(v_a(1), v_b(1), v_b(2), &
mg%lsf_boundary_value, dd, val)
box%fc(j, i, 2, i_fc) = grad_sign * fac * (val - v_a(2)) / &
box%fc(j, i, 2, i_fc) = grad_sign * fac * (v_a(2) - val) / &
(box%dr(2) * dd)
end do
end do
......@@ -2220,7 +2228,7 @@ contains
call lsf_dist_val(v_a(1), v_b(1), v_b(2), &
mg%lsf_boundary_value, dd, val)
box%fc(IJK, 1, i_fc) = grad_sign * fac * (val - v_a(2)) / &
box%fc(IJK, 1, i_fc) = grad_sign * fac * (v_a(2) - val) / &
(box%dr(1) * dd)
if (box%cc(j, i, k, i_lsf) > 0) then
......@@ -2235,7 +2243,7 @@ contains
call lsf_dist_val(v_a(1), v_b(1), v_b(2), &
mg%lsf_boundary_value, dd, val)
box%fc(j, i, k, 2, i_fc) = grad_sign * fac * (val - v_a(2)) / &
box%fc(j, i, k, 2, i_fc) = grad_sign * fac * (v_a(2) - val) / &
(box%dr(2) * dd)
if (box%cc(j, k, i, i_lsf) > 0) then
......@@ -2250,7 +2258,7 @@ contains
call lsf_dist_val(v_a(1), v_b(1), v_b(2), &
mg%lsf_boundary_value, dd, val)
box%fc(j, k, i, 3, i_fc) = grad_sign * fac * (val - v_a(2)) / &
box%fc(j, k, i, 3, i_fc) = grad_sign * fac * (v_a(2) - val) / &
(box%dr(3) * dd)
end do
end do
......@@ -2279,7 +2287,8 @@ contains
end do
if (max_lsf * min_lsf >= 0) then
print *, "Make the coarse grid finer?"
print *, "Min/max level set function:", min_lsf, max_lsf
print *, "Perhaps make the coarse grid finer?"
error stop "Level set function does not change sign on coarse grid"
end if
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment