Commit 003ab6ed authored by Jannis Teunissen's avatar Jannis Teunissen

Multigrid: add level-set function for internal boundary conditions

parent 2463b01f
......@@ -55,3 +55,4 @@
documentation/html/*
external_libraries/*
tests/results/*
/examples/electrode_example_?d
......@@ -10,7 +10,8 @@ include $(AF_DIR)/src/makerules.make
PROGS_XD := random_refinement poisson_basic poisson_benchmark advection \
computational_domain boundary_conditions poisson_neumann particles_to_grid \
ghostcell_benchmark particles_gravity poisson_coarse_solver implicit_diffusion \
poisson_dielectric poisson_helmholtz reaction_diffusion dielectric_surface advection_v2
poisson_dielectric poisson_helmholtz reaction_diffusion dielectric_surface \
advection_v2 electrode_example
PROGS_1D := $(PROGS_XD:%=%_1d)
......
#include "../src/cpp_macros.h"
!> \example dielectric_surface.f90
!>
!> Example showing how to include a dielectric surface
program dielectric_surface
use m_af_all
use m_dielectric
implicit none
integer, parameter :: box_size = 8
integer, parameter :: n_iterations = 10
integer :: i_phi
integer :: i_rhs
integer :: i_tmp
integer :: i_lsf
type(af_t) :: tree
type(ref_info_t) :: ref_info
integer :: n, mg_iter, coord, n_args
real(dp) :: residu
character(len=100) :: fname
type(mg_t) :: mg
call af_add_cc_variable(tree, "phi", ix=i_phi)
call af_add_cc_variable(tree, "rhs", ix=i_rhs)
call af_add_cc_variable(tree, "tmp", ix=i_tmp)
call af_add_cc_variable(tree, "lsf", ix=i_lsf)
call af_set_cc_methods(tree, i_lsf, af_bc_neumann_zero)
call af_set_cc_methods(tree, i_rhs, af_bc_neumann_zero)
! If an argument is given, switch to cylindrical coordinates in 2D
n_args = command_argument_count()
if (NDIM == 2 .and. n_args == 1) then
coord = af_cyl
else
coord = af_xyz
end if
! Initialize tree
call af_init(tree, & ! Tree to initialize
box_size, & ! A box contains box_size**DIM cells
[DTIMES(1.0_dp)], &
2 * [DTIMES(box_size)], &
coord=coord)
do n = 1, 4
call af_adjust_refinement(tree, ref_routine, ref_info)
if (ref_info%n_add == 0) exit
end do
call af_loop_box(tree, set_init_cond)
mg%i_phi = i_phi
mg%i_rhs = i_rhs
mg%i_tmp = i_tmp
mg%i_lsf = i_lsf
mg%sides_bc => af_bc_dirichlet_zero
mg%lsf_boundary_value = 1.0_dp
call mg_init(tree, mg)
do mg_iter = 1, n_iterations
call mg_fas_fmg(tree, mg, .true., mg_iter>1)
! Determine the minimum and maximum residual and error
call af_tree_maxabs_cc(tree, i_tmp, residu)
write(*, "(I8,Es14.5)") mg_iter, residu
write(fname, "(A,I0)") "electrode_example_" // DIMNAME // "_", mg_iter
call af_write_silo(tree, trim(fname), dir="output")
end do
contains
! Set refinement flags for box
subroutine ref_routine(box, cell_flags)
type(box_t), intent(in) :: box
integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
if (box%lvl < 5 .and. box%r_min(1) < 0.5_dp) then
cell_flags(DTIMES(:)) = af_do_ref
else
cell_flags(DTIMES(:)) = af_keep_ref
end if
end subroutine ref_routine
! This routine sets the initial conditions for each box
subroutine set_init_cond(box)
type(box_t), intent(inout) :: box
integer :: IJK, nc
real(dp) :: rr(NDIM), r0(NDIM), r1(NDIM)
real(dp) :: radius
nc = box%n_cell
! Start and end point of line
r0(:) = 0.4_dp
r1(:) = 0.6_dp
radius = 0.02_dp
do KJI_DO(0,nc+1)
rr = af_r_cc(box, [IJK])
box%cc(IJK, i_lsf) = dist_vec_line(rr, r0, r1, NDIM) - radius
end do; CLOSE_DO
end subroutine set_init_cond
!> Compute distance vector between point and its projection onto a line
!> between r0 and r1
function dist_vec_line(r, r0, r1, n_dim) result(dist)
integer, intent(in) :: n_dim
real(dp), intent(in) :: r(n_dim), r0(n_dim), r1(n_dim)
real(dp) :: dist_vec(n_dim), frac
real(dp) :: dist, line_len2
line_len2 = sum((r1 - r0)**2)
frac = sum((r - r0) * (r1 - r0))
if (frac <= 0.0_dp) then
dist_vec = r - r0
else if (frac >= line_len2) then
dist_vec = r - r1
else
dist_vec = r - (r0 + frac/line_len2 * (r1 - r0))
end if
dist = norm2(dist_vec)
end function dist_vec_line
end program dielectric_surface
This diff is collapsed.
......@@ -91,7 +91,7 @@ contains
call hypre_create_vector(mg%csolver%grid, nx, mg%csolver%phi)
call hypre_create_vector(mg%csolver%grid, nx, mg%csolver%rhs)
if (tree%coord_t == af_cyl) then
if (tree%coord_t == af_cyl .or. mg%i_lsf /= -1) then
! The symmetry option does not seem to work well with axisymmetric problems
mg%csolver%symmetric = 0
end if
......@@ -123,12 +123,17 @@ contains
stencil_size, stencil_offsets(:, stencil_ix), mg%csolver%symmetric)
allocate(mg%csolver%bc_to_rhs(nc**(NDIM-1), af_num_neighbors, n_boxes))
allocate(mg%csolver%lsf_fac(DTIMES(nc), n_boxes))
allocate(coeffs(stencil_size, tree%n_cell**NDIM))
mg%csolver%bc_to_rhs = 0.0_dp
mg%csolver%lsf_fac = 0.0_dp
coeffs = 0.0_dp
do n = 1, size(tree%lvls(1)%ids)
id = tree%lvls(1)%ids(n)
call mg%box_stencil(tree%boxes(id), mg, full_coeffs, &
mg%csolver%bc_to_rhs(:, :, n))
mg%csolver%bc_to_rhs(:, :, n), mg%csolver%lsf_fac(DTIMES(:), n))
cnt = 0
do KJI_DO(1, nc)
......@@ -263,6 +268,11 @@ contains
end if
end do
! Add contribution of level-set function to rhs
if (mg%i_lsf /= -1) then
tmp = tmp + mg%csolver%lsf_fac(DTIMES(:), n) * mg%lsf_boundary_value
end if
ilo = (tree%boxes(id)%ix - 1) * nc + 1
ihi = ilo + nc - 1
call HYPRE_StructVectorSetBoxValues(mg%csolver%rhs, ilo, ihi, &
......
......@@ -29,6 +29,9 @@ module m_mg_types
!> Stores coefficient to convert boundary conditions to the right-hand side
real(dp), allocatable :: bc_to_rhs(:, :, :)
!> Stores coefficients to use with level set function
real(dp), allocatable :: lsf_fac(DTIMES(:), :)
integer :: symmetric = 1
integer :: solver_type = -1
integer :: max_iterations = 50
......@@ -45,7 +48,6 @@ module m_mg_types
integer :: i_eps = -1 !< Optional variable (diel. permittivity)
integer :: i_lsf = -1 !< Optional variable for level set function
integer :: i_bval = -1 !< Optional variable for boundary value
integer :: n_cycle_down = -1 !< Number of relaxation cycles in downward sweep
integer :: n_cycle_up = -1 !< Number of relaxation cycles in upward sweep
......@@ -57,6 +59,9 @@ module m_mg_types
!> Store lambda^2 for Helmholtz equations (L phi - lamda phi = f)
real(dp) :: helmholtz_lambda = 0.0_dp
!> Boundary value for level set function
real(dp) :: lsf_boundary_value = 0.0_dp
!> Routine to call for filling ghost cells near physical boundaries
procedure(af_subr_bc), pointer, nopass :: sides_bc => null()
......@@ -113,12 +118,13 @@ module m_mg_types
type(mg_t), intent(in) :: mg !< Multigrid options
end subroutine mg_box_rstr
subroutine mg_box_stencil(box, mg, stencil, bc_to_rhs)
subroutine mg_box_stencil(box, mg, stencil, bc_to_rhs, lsf_fac)
import
type(box_t), intent(in) :: box
type(mg_t), intent(in) :: mg
real(dp), intent(inout) :: stencil(2*NDIM+1, DTIMES(box%n_cell))
real(dp), intent(inout) :: bc_to_rhs(box%n_cell**(NDIM-1), af_num_neighbors)
real(dp), intent(inout) :: lsf_fac(DTIMES(box%n_cell))
end subroutine mg_box_stencil
end interface
......
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