Commit 9d80af8e authored by Jannis Teunissen's avatar Jannis Teunissen

Major update: add 1D version of Afivo framework

Some of the updates:
- All library functionality has been included in 1D (e.g. multigrid, particles)
- Silo output now uses curve format (standard multimesh did not work in Visit)
- Hypre uses PCG method because PFMG/SMG work only in 2D/3D
- Examples have been extended to 1D
parent df9d63c5
.DS_Store
*.a
*.dat
*.mod
*.o
*.vtu
*.silo
*.a
src/m_a2_*
src/m_a3_*
/examples/poisson_basic_2d*
/examples/poisson_basic_3d*
/examples/advection_2d*
/examples/advection_3d*
/examples/boundary_conditions_2d*
/examples/boundary_conditions_3d*
/examples/computational_domain_2d*
/examples/computational_domain_3d*
/examples/implicit_diffusion_2d*
/examples/implicit_diffusion_3d*
/examples/poisson_benchmark_2d*
/examples/poisson_benchmark_3d*
/examples/drift_diffusion_2d*
/examples/drift_diffusion_3d*
/examples/random_refinement_2d*
/examples/random_refinement_3d*
/examples/poisson_neumann_2d*
/examples/poisson_neumann_3d*
*.vtk
*.vtu
./examples/*.curve
./examples/*.out
.DS_Store
/documentation/auto/
/documentation/doc/
/examples/advection_?d*
/examples/advection_v2_?d
/examples/amr_solid_body_rotation
/examples/boundary_conditions_?d*
/examples/casper_test
/examples/computational_domain_?d*
/examples/dielectric_?d
/examples/dielectric_?d.f90
/examples/dielectric_surface_?d
/examples/drift_diffusion_?d*
/examples/euler_gas_dynamics
/examples/ghostcell_benchmark_?d
/examples/ghostcell_benchmark_?d.f90
/examples/implicit_diffusion_?d*
/examples/particles_gravity_?d
/examples/particles_to_grid_?d
/examples/particles_to_grid_?d.f90
/examples/poisson_basic_?d*
/examples/poisson_benchmark_?d*
/examples/poisson_coarse_solver_?d
/examples/poisson_cyl*
/examples/poisson_cyl_dielectric*
/tests/test_morton
/tests/test_reduction_2d
/documentation/doc/
/examples/poisson_dielectric_?d
/examples/poisson_div_cleaning
/examples/poisson_helmholtz_?d
/examples/poisson_helmholtz_?d.f90
/examples/poisson_helmholtz_cyl
/examples/poisson_neumann_?d*
/examples/random_refinement_?d*
/examples/reaction_diffusion_?d
/examples/simple_streamer
/examples/simple_streamer_?d
/examples/solid_body_rotation
/tests/test_ghostcell
/tests/test_ghostcell_?d
/tests/test_init
/tests/test_morton
/tests/test_morton_?d
/tests/test_reduction_?d
/tests/test_refinement
/tests/test_types_2d_3d
tests/results/*
/examples/simple_streamer_2d
external_libraries/*
/tests/test_refinement_?d
/tests/test_types_?d
documentation/html/*
*.dat
*.vtk
/examples/particles_to_grid_2d
/examples/particles_to_grid_3d
/examples/particles_to_grid_2d.f90
/examples/particles_to_grid_3d.f90
/examples/casper_test
/examples/dielectric_2d
/examples/dielectric_2d.f90
/examples/dielectric_3d
/examples/dielectric_3d.f90
/examples/ghostcell_benchmark_2d
/examples/ghostcell_benchmark_2d.f90
/examples/ghostcell_benchmark_3d
/examples/ghostcell_benchmark_3d.f90
/examples/poisson_helmholtz_2d
/examples/poisson_helmholtz_2d.f90
/examples/poisson_helmholtz_3d
/examples/poisson_helmholtz_3d.f90
/examples/particles_gravity_2d
/examples/particles_gravity_3d
/examples/poisson_dielectric_2d
/examples/poisson_dielectric_3d
/examples/poisson_helmholtz_cyl
/tests/test_ghostcell_2d
/tests/test_ghostcell_3d
/tests/test_morton_2d
/tests/test_morton_3d
/tests/test_reduction_3d
/tests/test_refinement_2d
/tests/test_refinement_3d
/tests/test_types_2d
/tests/test_types_3d
/examples/simple_streamer
/examples/poisson_coarse_solver_2d
/examples/poisson_coarse_solver_3d
/examples/poisson_div_cleaning
/examples/reaction_diffusion_2d
/examples/reaction_diffusion_3d
/examples/dielectric_surface_2d
/examples/dielectric_surface_3d
/documentation/auto/
./examples/*.curve
./examples/*.out
/examples/euler_gas_dynamics
/examples/solid_body_rotation
/examples/amr_solid_body_rotation
/examples/advection_v2_2d
/examples/advection_v2_3d
external_libraries/*
tests/results/*
SRC_DIRS := lib_2d lib_3d examples tests
SRC_DIRS := lib_1d lib_2d lib_3d examples tests
# Directories with altered names (useful for cleaning)
CLEANSRC := $(SRC_DIRS:%=clean-%)
......@@ -30,4 +30,4 @@ external_libraries/hypre:
# Dependency information
$(SRC_DIRS): external_libraries/silo external_libraries/hypre
examples tests: lib_2d lib_3d
examples tests: lib_1d lib_2d lib_3d
......@@ -12,12 +12,14 @@ PROGS_XD := random_refinement poisson_basic poisson_benchmark advection \
ghostcell_benchmark particles_gravity poisson_coarse_solver implicit_diffusion \
poisson_dielectric poisson_helmholtz reaction_diffusion dielectric_surface advection_v2
PROGS_1D := $(PROGS_XD:%=%_1d)
PROGS_2D := $(PROGS_XD:%=%_2d) poisson_cyl poisson_cyl_dielectric \
simple_streamer poisson_cyl_analytic poisson_helmholtz_cyl solid_body_rotation amr_solid_body_rotation euler_gas_dynamics
PROGS_3D := $(PROGS_XD:%=%_3d) poisson_div_cleaning
PROGS := $(PROGS_2D) $(PROGS_3D)
PROGS := $(PROGS_1D) $(PROGS_2D) $(PROGS_3D)
all: $(PROGS)
......@@ -42,7 +44,12 @@ run_3d:
$(OUTDIR):
mkdir -p $@
# Set flags for 2d and 3d version of programs
# Set flags
$(PROGS_1D): $(AF_DIR)/lib_1d/libafivo.a | $(OUTDIR)
$(PROGS_1D): INCDIRS+=$(AF_DIR)/lib_1d
$(PROGS_1D): LIBDIRS+=$(AF_DIR)/lib_1d
$(PROGS_1D): FFLAGS += -DNDIM=1
$(PROGS_2D): $(AF_DIR)/lib_2d/libafivo.a | $(OUTDIR)
$(PROGS_2D): INCDIRS+=$(AF_DIR)/lib_2d
$(PROGS_2D): LIBDIRS+=$(AF_DIR)/lib_2d
......@@ -54,6 +61,9 @@ $(PROGS_3D): LIBDIRS+=$(AF_DIR)/lib_3d
$(PROGS_3D): FFLAGS += -DNDIM=3
# How to create executables
%_1d: %.f90
$(FC) -o $@ $^ $(FFLAGS) $(addprefix -I,$(INCDIRS)) \
$(addprefix -L,$(LIBDIRS)) $(addprefix -l,$(LIBS))
%_2d: %.f90
$(FC) -o $@ $^ $(FFLAGS) $(addprefix -I,$(INCDIRS)) \
$(addprefix -L,$(LIBDIRS)) $(addprefix -l,$(LIBS))
......
......@@ -13,7 +13,7 @@ program advection
integer :: i_phi_old
integer :: i_err
integer :: i_flux
integer, parameter :: sol_type = 1
integer, parameter :: sol_type = 2
! Set coord_type to af_cyl to test conservation in cyl. coordinates
integer, parameter :: coord_type = af_xyz
real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
......@@ -55,9 +55,8 @@ program advection
dt_adapt = 0.01_dp
dt_output = 0.5_dp
end_time = 5.0_dp
velocity(:) = 0.0_dp
velocity(:) = -0.5_dp
velocity(1) = 1.0_dp
velocity(2) = -1.0_dp
! Set up the initial conditions
call system_clock(t_start,count_rate)
......@@ -189,7 +188,10 @@ contains
nc = box%n_cell
do KJI_DO(1,nc)
#if NDIM == 2
#if NDIM == 1
diff = abs(box%dr(1) * (box%cc(i+1, i_phi) + &
box%cc(i-1, i_phi) - 2 * box%cc(i, i_phi)))
#elif NDIM == 2
diff = abs(box%dr(1) * (box%cc(i+1, j, i_phi) + &
box%cc(i-1, j, i_phi) - 2 * box%cc(i, j, i_phi)) + &
box%dr(2) * (box%cc(i, j+1, i_phi) + &
......@@ -251,7 +253,11 @@ contains
select case (sol_type)
case (1)
#if NDIM > 1
sol = sin(0.5_dp * rr_t(1))**8 * cos(0.5_dp * rr_t(2))**8
#else
sol = sin(0.5_dp * rr_t(1))**8
#endif
case (2)
rr_t = modulo(rr_t, domain_len) / domain_len
if (norm2(rr_t - 0.5_dp) < 0.1_dp) then
......@@ -278,7 +284,11 @@ contains
call af_gc2_box(tree, id, [i_phi], cc)
#if NDIM == 2
#if NDIM == 1
v(:, 1) = velocity(1)
call flux_koren_1d(cc(DTIMES(:), 1), v, nc, 2)
tree%boxes(id)%fc(:, :, i_phi) = v
#elif NDIM == 2
v(:, :, 1) = velocity(1)
v(:, :, 2) = velocity(2)
......@@ -308,7 +318,13 @@ contains
nc = box%n_cell
inv_dr = 1/box%dr
#if NDIM == 2
#if NDIM == 1
forall (i = 1:nc)
box%cc(i, i_phi) = box%cc(i, i_phi) + dt(1) * ( &
inv_dr(1) * &
(box%fc(i, 1, i_phi) - box%fc(i+1, 1, i_phi)))
end forall
#elif NDIM == 2
if (coord_type == af_cyl) then
call af_cyl_flux_factors(box, rfac)
do j = 1, nc
......
......@@ -11,7 +11,7 @@ program advection
integer :: i_phi
integer :: i_err
integer :: i_flux
integer, parameter :: sol_type = 1
integer, parameter :: sol_type = 2
real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
type(af_t) :: tree
......@@ -48,9 +48,8 @@ program advection
dt_amr = 0.01_dp
dt_output = 0.5_dp
end_time = 5.0_dp
velocity(:) = 0.0_dp
velocity(:) = -0.5_dp
velocity(1) = 1.0_dp
velocity(2) = -1.0_dp
! Set up the initial conditions
refine_steps=0
......@@ -134,7 +133,10 @@ contains
nc = box%n_cell
do KJI_DO(1,nc)
#if NDIM == 2
#if NDIM == 1
diff = abs(box%dr(1) * (box%cc(i+1, i_phi) + &
box%cc(i-1, i_phi) - 2 * box%cc(i, i_phi)))
#elif NDIM == 2
diff = abs(box%dr(1) * (box%cc(i+1, j, i_phi) + &
box%cc(i-1, j, i_phi) - 2 * box%cc(i, j, i_phi)) + &
box%dr(2) * (box%cc(i, j+1, i_phi) + &
......@@ -196,7 +198,11 @@ contains
select case (sol_type)
case (1)
#if NDIM > 1
sol = sin(0.5_dp * rr_t(1))**8 * cos(0.5_dp * rr_t(2))**8
#else
sol = sin(0.5_dp * rr_t(1))**8
#endif
case (2)
rr_t = modulo(rr_t, domain_len) / domain_len
if (norm2(rr_t - 0.5_dp) < 0.1_dp) then
......
......@@ -41,7 +41,7 @@ program boundary_conditions
call af_gc_tree(tree, [i_phi])
write(fname, "(A,I0)") "boundary_conditions_" // DIMNAME // "_", iter
call af_write_vtk(tree, trim(fname), dir="output")
call af_write_silo(tree, trim(fname), dir="output")
end do
contains
......@@ -61,7 +61,11 @@ contains
nc = box%n_cell
do KJI_DO(1,nc)
#if NDIM == 2
#if NDIM == 1
tmp(i) = 0.5_dp * ( &
box%cc(i+1, i_phi) + &
box%cc(i-1, i_phi))
#elif NDIM == 2
tmp(i, j) = 0.25_dp * ( &
box%cc(i+1, j, i_phi) + &
box%cc(i-1, j, i_phi) + &
......@@ -78,7 +82,9 @@ contains
#endif
end do; CLOSE_DO
box%cc(DTIMES(1:nc), i_phi) = tmp(DTIMES(:))
! Average new and old value
box%cc(DTIMES(1:nc), i_phi) = 0.5_dp * (&
tmp(DTIMES(:)) + box%cc(DTIMES(1:nc), i_phi))
end subroutine average_phi
!> [boundary_method]
......@@ -95,39 +101,12 @@ contains
! Below the solution is specified in the approriate ghost cells
select case (nb)
#if NDIM == 2
case (af_neighb_lowx) ! Lower-x direction
bc_type = af_bc_dirichlet
bc_val = 1.0_dp
case (af_neighb_highx) ! Higher-x direction
case default
bc_type = af_bc_neumann
bc_val = 0.0_dp
case (af_neighb_lowy) ! Lower-y direction
bc_type = af_bc_dirichlet
bc_val = 1.0_dp
case (af_neighb_highy) ! Higher-y direction
bc_type = af_bc_neumann
bc_val = 0.0_dp
#elif NDIM == 3
case (af_neighb_lowx) ! Lower-x direction
bc_type = af_bc_dirichlet
bc_val = 1.0_dp
case (af_neighb_highx) ! Higher-x direction
bc_type = af_bc_neumann
bc_val = 0.0_dp
case (af_neighb_lowy) ! Lower-y direction
bc_type = af_bc_dirichlet
bc_val = 1.0_dp
case (af_neighb_highy) ! Higher-y direction
bc_type = af_bc_neumann
bc_val = 0.0_dp
case (af_neighb_lowz) ! Lower-z direction
bc_type = af_bc_neumann
bc_val = 0.0_dp
case (af_neighb_highz) ! Higher-z direction
bc_type = af_bc_neumann
bc_val = 0.0_dp
#endif
end select
end subroutine boundary_method
......
......@@ -18,18 +18,20 @@ program computational_domain
domain_size = 1.0_dp * grid_size
call af_add_cc_variable(tree, "phi")
call af_init(tree, n_cell, domain_size, grid_size)
call af_write_vtk(tree, "computational_domain_" // DIMNAME // "_1", dir="output")
call af_init(tree, n_cell, domain_size, grid_size, mem_limit_gb=0.1_dp)
call af_write_silo(tree, "computational_domain_" // DIMNAME // "_1", dir="output")
call af_destroy(tree)
! Create mesh 2: two boxes along y-direction
grid_size(:) = n_cell
#if NDIM > 1
grid_size(2) = 2 * n_cell
#endif
domain_size = 1.0_dp * grid_size
call af_add_cc_variable(tree, "phi")
call af_init(tree, n_cell, domain_size, grid_size)
call af_write_vtk(tree, "computational_domain_" // DIMNAME // "_2", dir="output")
call af_init(tree, n_cell, domain_size, grid_size, mem_limit_gb=0.1_dp)
call af_write_silo(tree, "computational_domain_" // DIMNAME // "_2", dir="output")
call af_destroy(tree)
! Create mesh 3: Two boxes along x-direction that are fully periodic
......@@ -39,8 +41,8 @@ program computational_domain
domain_size = 1.0_dp * grid_size
call af_add_cc_variable(tree, "phi")
call af_init(tree, n_cell, domain_size, grid_size, periodic=periodic)
call af_write_vtk(tree, "computational_domain_" // DIMNAME // "_3", dir="output")
call af_init(tree, n_cell, domain_size, grid_size, mem_limit_gb=0.1_dp)
call af_write_silo(tree, "computational_domain_" // DIMNAME // "_3", dir="output")
call af_destroy(tree)
end program computational_domain
......@@ -43,7 +43,9 @@ program dielectric_surface
! Include both cell-centered and face-centered electric field for testing
call af_add_cc_variable(tree, "fld_cc_x", ix=i_fld_cc(1))
#if NDIM > 1
call af_add_cc_variable(tree, "fld_cc_y", ix=i_fld_cc(2))
#endif
#if NDIM == 3
call af_add_cc_variable(tree, "fld_cc_z", ix=i_fld_cc(3))
#endif
......
......@@ -16,7 +16,9 @@ program implicit_diffusion
real(dp), parameter :: domain_len = 2 * acos(-1.0_dp)
real(dp), parameter :: diffusion_coeff = 1.0_dp
#if NDIM == 2
#if NDIM == 1
real(dp), parameter :: solution_modes(NDIM) = [1]
#elif NDIM == 2
real(dp), parameter :: solution_modes(NDIM) = [1, 1]
#elif NDIM == 3
real(dp), parameter :: solution_modes(NDIM) = [1, 1, 0]
......
......@@ -47,7 +47,9 @@ program particles_gravity
call af_add_cc_variable(tree, "rho", ix=i_rho)
call af_add_cc_variable(tree, "tmp", ix=i_tmp)
call af_add_cc_variable(tree, "fx", ix=i_f(1))
#if NDIM > 1
call af_add_cc_variable(tree, "fy", ix=i_f(2))
#endif
#if NDIM == 3
call af_add_cc_variable(tree, "fz", ix=i_f(3))
#endif
......@@ -241,7 +243,10 @@ contains
fac = -0.5_dp / box%dr
do KJI_DO(1, nc)
#if NDIM == 2
#if NDIM == 1
box%cc(i, i_f(1)) = fac(1) * &
(box%cc(i+1, i_phi) - box%cc(i-1, i_phi))
#elif NDIM == 2
box%cc(i, j, i_f(1)) = fac(1) * &
(box%cc(i+1, j, i_phi) - box%cc(i-1, j, i_phi))
box%cc(i, j, i_f(2)) = fac(2) * &
......
......@@ -9,11 +9,12 @@ program poisson_basic
implicit none
integer, parameter :: box_size = 8
integer, parameter :: box_size = 64
integer, parameter :: n_iterations = 10
integer :: domain_size(NDIM)
real(dp) :: domain_len(NDIM)
integer :: i_phi
integer :: i_sol
integer :: i_rhs
integer :: i_err
integer :: i_tmp
......@@ -43,6 +44,7 @@ program poisson_basic
!> [af_init]
call af_add_cc_variable(tree, "phi", ix=i_phi)
call af_add_cc_variable(tree, "sol", ix=i_sol)
call af_add_cc_variable(tree, "rhs", ix=i_rhs)
call af_add_cc_variable(tree, "err", ix=i_err)
call af_add_cc_variable(tree, "tmp", ix=i_tmp)
......@@ -180,6 +182,7 @@ contains
box%cc(IJK, i_rhs) = gauss_laplacian(gs, rr)
call gauss_gradient(gs, rr, grad)
box%cc(IJK, i_gradx) = grad(1)
box%cc(IJK, i_sol) = gauss_value(gs, rr)
end do; CLOSE_DO
end subroutine set_initial_condition
!> [set_initial_condition]
......@@ -196,7 +199,10 @@ contains
do KJI_DO(1,nc)
rr = af_r_cc(box, [IJK])
box%cc(IJK, i_err) = box%cc(IJK, i_phi) - gauss_value(gs, rr)
#if NDIM == 2
#if NDIM == 1
gradx = 0.5_dp * (box%cc(i+1, i_phi) - &
box%cc(i-1, i_phi)) / box%dr(1)
#elif NDIM == 2
gradx = 0.5_dp * (box%cc(i+1, j, i_phi) - &
box%cc(i-1, j, i_phi)) / box%dr(1)
#elif NDIM == 3
......
......@@ -187,7 +187,10 @@ contains
do KJI_DO(1,nc)
rr = af_r_cc(box, [IJK])
box%cc(IJK, i_err) = box%cc(IJK, i_phi) - gauss_value(gs, rr)
#if NDIM == 2
#if NDIM == 1
gradx = 0.5_dp * (box%cc(i+1, i_phi) - &
box%cc(i-1, i_phi)) / box%dr(1)
#elif NDIM == 2
gradx = 0.5_dp * (box%cc(i+1, j, i_phi) - &
box%cc(i-1, j, i_phi)) / box%dr(1)
#elif NDIM == 3
......
......@@ -140,7 +140,11 @@ contains
rr = af_r_cc(box, [IJK])
! Set the values at each cell according to some function
box%cc(IJK, i_phi) = sin(0.5_dp * rr(1))**2 * cos(rr(2))**2
#if NDIM > 1
box%cc(IJK, i_phi) = sin(0.5_dp * rr(1))**2 * cos(rr(2))**2
#else
box%cc(IJK, i_phi) = sin(0.5_dp * rr(1))**2
#endif
end do; CLOSE_DO
end subroutine set_init_cond
......
......@@ -397,7 +397,9 @@ contains
associate (cc => box%cc, n => i_in)
do KJI_DO(1, nc)
#if NDIM == 2
#if NDIM == 1
lpl(i) = idr2(1) * (cc(i-1, n) + cc(i+1, n) - 2 * cc(i, n))
#elif NDIM == 2
lpl(i, j) = &
idr2(1) * (cc(i-1, j, n) + cc(i+1, j, n) - 2 * cc(i, j, n)) + &
idr2(2) * (cc(i, j-1, n) + cc(i, j+1, n) - 2 * cc(i, j, n))
......
NDIM := 1
LIB := libafivo.a
.PHONY: all clean
all: $(LIB)
vpath %.f90 ../src
include ../src/definitions.make
include ../src/makerules.make
$(OBJS): FFLAGS += -DNDIM=$(NDIM)
$(LIB): $(OBJS)
$(RM) $@
$(AR) rcs $@ $^
clean:
$(RM) $(LIB) *.o *.mod
......@@ -8,7 +8,7 @@ vpath %.f90 ../src
include ../src/definitions.make
include ../src/makerules.make
$(OBJS): FFLAGS += -DNDIM=2
$(OBJS): FFLAGS += -DNDIM=$(NDIM)
$(LIB): $(OBJS)
$(RM) $@
......
NDIM := 3
LIB := libafivo.a
.PHONY: all clean
......@@ -8,7 +9,7 @@ vpath %.f90 ../src
include ../src/definitions.make
include ../src/makerules.make
$(OBJS): FFLAGS += -DNDIM=3
$(OBJS): FFLAGS += -DNDIM=$(NDIM)
$(LIB): $(OBJS)
$(RM) $@
......
#if NDIM == 2
#ifdef __GFORTRAN__
#define PASTE(a) a
#define CONCAT(a,b) PASTE(a)b
#else
#define PASTE(a) a ## b
#define CONCAT(a,b) PASTE(a,b)