Test of the level-set functionality of the Poisson solver.
Test of the level-set functionality of the Poisson solver In this example, multiple irregular boundaries can be defined, for example a sphere or a heart shape. For the case of a sphere, an analytic solution is provided in 1d, 2d and 3d.
1#include "../src/cpp_macros.h"
2
3
4
5
6
7
8
9
10
11
12
13
14program poisson_lsf_test
17 use omp_lib
18
19 implicit none
20
21 integer :: box_size = 8
22 integer :: n_iterations = 10
23 integer :: max_refine_level = 4
24 integer :: min_refine_level = 2
25 integer :: i_phi
26 integer :: i_rhs
27 integer :: i_tmp
28 integer :: i_lsf
29 integer :: i_lsf_mask
30 integer :: i_error
31 integer :: i_field
32 integer :: i_field_norm
33
34 logical :: cylindrical = .false.
35 logical :: write_numpy = .false.
36 integer :: refinement_type = 1
37
38
39 integer :: shape = 1
40
41 real(dp) :: sharpness_t = 4.0_dp
42 real(dp) :: boundary_value = 1.0_dp
43 real(dp) :: solution_coeff = 1.0_dp
44 real(dp) :: solution_radius = 0.25_dp
45 real(dp) :: solution_r0(NDIM) = 0.5_dp
46 real(dp) :: solution_smallest_width = 1e-3_dp
47
48 type(af_t) :: tree
49 type(ref_info_t) :: ref_info
50 integer :: n, mg_iter, coord
51 real(dp) :: residu, max_error, max_field
52 character(len=150) :: fname
53 character(len=20) :: output_suffix = ""
54 character(len=20) :: prolongation = "auto"
55 type(mg_t) :: mg
56 type(CFG_t) :: cfg
57 real(dp) :: error_squared, rmse
58 real(dp) :: mem_limit_gb = 8.0_dp
59 real(dp) :: t0, t_sum
60 logical :: write_output = .true.
61
62 call af_add_cc_variable(tree, "phi", ix=i_phi)
63 call af_add_cc_variable(tree, "rhs", ix=i_rhs)
64 call af_add_cc_variable(tree, "tmp", ix=i_tmp)
65 call af_add_cc_variable(tree, "lsf", ix=i_lsf)
66 call af_add_cc_variable(tree, "lsf_mask", ix=i_lsf_mask)
67 call af_add_cc_variable(tree, "error", ix=i_error)
68 call af_add_cc_variable(tree, "field_norm", ix=i_field_norm)
69 call af_add_fc_variable(tree, "field", ix=i_field)
70
71 call af_set_cc_methods(tree, i_lsf, funcval=set_lsf)
72 call af_set_cc_methods(tree, i_lsf_mask, af_bc_neumann_zero)
73 call af_set_cc_methods(tree, i_field_norm, af_bc_neumann_zero)
74
75 call cfg_update_from_arguments(cfg)
76
78 "Whether to use cylindrical coordinates")
79 call cfg_add_get(cfg,
"write_output", write_output, &
80 "Whether to write Silo output")
82 "Whether to write Numpy output")
83 call cfg_add_get(cfg,
"output_suffix", output_suffix, &
84 "Suffix to add to output")
85 call cfg_add_get(cfg,
"max_refine_level", max_refine_level, &
86 "Maximum refinement level")
87 call cfg_add_get(cfg,
"min_refine_level", min_refine_level, &
88 "Minimum refinement level")
89 call cfg_add_get(cfg,
"refinement_type", refinement_type, &
90 "Type of refinement (1: uniform, 2: only boundary)")
92 "Index of the shape used")
94 "Size of grid boxes")
95 call cfg_add_get(cfg,
"n_iterations", n_iterations, &
96 "Number of multigrid iterations")
97 call cfg_add_get(cfg,
"n_cycle_up", mg%n_cycle_up, &
98 "Number of relaxation steps going up")
99 call cfg_add_get(cfg,
"n_cycle_down", mg%n_cycle_down, &
100 "Number of relaxation steps going down")
101 call cfg_add_get(cfg,
"solution%sharpness", sharpness_t, &
102 "Sharpness of solution (for some cases)")
103 call cfg_add_get(cfg,
"solution%r0", solution_r0, &
104 "Origin of solution")
105 call cfg_add_get(cfg,
"solution%radius", solution_radius, &
106 "Solution radius")
107 call cfg_add_get(cfg,
"solution%coeff", solution_coeff, &
108 "Linear coefficient for solution")
109 call cfg_add_get(cfg,
"solution%smallest_width", solution_smallest_width, &
110 "Smallest width to resolve in the solution")
111 call cfg_add_get(cfg,
"boundary_value", boundary_value, &
112 "Value for Dirichlet boundary condition")
113 call cfg_add_get(cfg,
"mem_limit_gb", mem_limit_gb, &
114 "Memory limit (GByte)")
115 call cfg_add_get(cfg,
"prolongation", prolongation, &
116 "Prolongation method (linear, sparse, auto, custom)")
117
118 call cfg_check(cfg)
119
120
121 if (cylindrical) then
122 coord = af_cyl
123
124 solution_r0(1) = 0.0_dp
125 else
126 coord = af_xyz
127 end if
128
129 tree%mg_i_lsf = i_lsf
130 mg%i_phi = i_phi
131 mg%i_rhs = i_rhs
132 mg%i_tmp = i_tmp
133 mg%sides_bc => bc_solution
134 mg%lsf_boundary_value = boundary_value
135 mg%lsf_dist => mg_lsf_dist_gss
136 mg%lsf => get_lsf
137 mg%lsf_length_scale = solution_smallest_width
138
139 select case (prolongation)
140 case ("linear")
141 mg%prolongation_type = mg_prolong_linear
142 case ("sparse")
143 mg%prolongation_type = mg_prolong_sparse
144 case ("auto")
145 mg%prolongation_type = mg_prolong_auto
146 case ("custom")
147 mg%lsf_use_custom_prolongation = .true.
148 case default
149 error stop "Unknown prolongation method"
150 end select
151
152
153 call af_init(tree, &
154 box_size, &
155 [dtimes(1.0_dp)], &
156 [dtimes(box_size)], &
157 coord=coord, &
158 mem_limit_gb=mem_limit_gb)
159
160 call af_refine_up_to_lvl(tree, min_refine_level)
161 call mg_init(tree, mg)
162
163 do n = 1, 100
164 call mg_set_operators_tree(tree, mg)
165 call af_adjust_refinement(tree, ref_routine, ref_info, ref_buffer=1)
166 if (ref_info%n_add == 0) exit
167 end do
168
169
170 call af_loop_box(tree, set_lsf_mask)
171 call af_gc_tree(tree, [i_lsf_mask])
172
173 call af_print_info(tree)
174
175 t_sum = 0
176 write(*, "(A,A4,4A14)") "# ", "iter", "residu", "max_error", "rmse", "max_field"
177
178 do mg_iter = 1, n_iterations
179 t0 = omp_get_wtime()
180 call mg_fas_fmg(tree, mg, .true., mg_iter>1)
181 t_sum = t_sum + omp_get_wtime() - t0
182
183 call mg_compute_phi_gradient(tree, mg, i_field, 1.0_dp, i_field_norm)
184 call af_gc_tree(tree, [i_field_norm])
185 call af_loop_box(tree, set_error)
186
187
188 call af_tree_maxabs_cc(tree, i_tmp, residu)
189 call af_tree_maxabs_cc(tree, i_error, max_error)
190 call af_tree_maxabs_cc(tree, i_field_norm, max_field)
191 call af_tree_sum_cc(tree, i_error, error_squared, 2)
192 rmse = sqrt(error_squared/af_total_volume(tree))
193 write(*, "(A,i4,4E14.5)") "# ", mg_iter, residu, max_error, rmse, max_field
194
195 write(fname, "(A,I0)") "output/poisson_lsf_test_" // dimname // &
196 trim(output_suffix) // "_", mg_iter
197 if (write_output) then
198 call af_write_silo(tree, trim(fname))
199 end if
200
201 if (write_numpy) then
202 call af_write_numpy(tree, trim(fname)//".npz", &
203 ixs_cc=[i_phi, i_lsf_mask])
204 end if
205 end do
206
207 write(*, "(A,E14.5)") " seconds_per_FMG_cycle", t_sum/n_iterations
208
209 call af_stencil_print_info(tree)
210
211contains
212
213
214 subroutine ref_routine(box, cell_flags)
215 type(box_t), intent(in) :: box
216 integer, intent(out) :: cell_flags(DTIMES(box%n_cell))
217 real(dp) :: tmp(DTIMES(box%n_cell)), dr_min
218 integer :: nc, IJK, ix_mask, n, cix(NDIM)
219
220 nc = box%n_cell
221 cell_flags(dtimes(:)) = af_keep_ref
222
223 select case (refinement_type)
224 case (1)
225
226 if (box%lvl < max_refine_level) then
227 cell_flags(dtimes(:)) = af_do_ref
228 end if
229
230 case (2)
231
232 ix_mask = af_stencil_index(box, mg_lsf_mask_key)
233 if (ix_mask /= af_stencil_none .and. box%lvl < max_refine_level) then
234 do n = 1, size(box%stencils(ix_mask)%sparse_ix, 2)
235 cix = box%stencils(ix_mask)%sparse_ix(:, n)
236 cell_flags(dindex(cix)) = af_do_ref
237 end do
238 end if
239
240 case (3)
241
242 if (norm2(box%r_min - solution_r0) < solution_radius .and. &
243 box%lvl < max_refine_level) then
244 cell_flags(dtimes(:)) = af_do_ref
245 end if
246 case (4)
247
248 dr_min = minval(af_lvl_dr(tree, max_refine_level))
249
250 do kji_do(1, nc)
251 tmp(ijk) = max(1.0_dp, norm2(af_r_cc(box, [ijk]) - solution_r0) &
252 / solution_radius)
253 end do; close_do
254
255 if (minval(box%dr) > minval(tmp) * dr_min .and. box%lvl < max_refine_level) then
256 cell_flags(dtimes(:)) = af_do_ref
257 end if
258 case default
259 error stop "Invalid refinement_type"
260 end select
261 end subroutine ref_routine
262
263 real(dp) function solution(r)
264 real(dp), intent(in) :: r(NDIM)
265 real(dp) :: distance, lsf
266
267 select case (shape)
268 case (1)
269
270 distance = norm2(r-solution_r0) / solution_radius
271
272
273 if (distance < 1.0_dp) then
274 solution = boundary_value
275 else if (ndim == 1) then
276 solution = boundary_value + solution_coeff * (distance - 1.0_dp)
277 else if (ndim == 2 .and. coord == af_xyz) then
278 solution = boundary_value + solution_coeff * log(distance)
279 else
280 solution = boundary_value + solution_coeff * (1 - 1/distance)
281 end if
282 case (5)
283
284 lsf = get_lsf(r)
285 if (lsf <= 0.0_dp) then
286 solution = boundary_value
287 else
288 solution = boundary_value * exp(-lsf)
289 end if
290 case default
291 solution = 0.0_dp
292 end select
293 end function solution
294
295
296 subroutine set_lsf(box, iv)
297 type(box_t), intent(inout) :: box
298 integer, intent(in) :: iv
299 integer :: IJK, nc
300 real(dp) :: rr(NDIM), norm_dr
301
302 nc = box%n_cell
303 norm_dr = norm2(box%dr)
304
305 do kji_do(0,nc+1)
306 rr = af_r_cc(box, [ijk])
307 box%cc(ijk, iv) = get_lsf(rr)
308 end do; close_do
309 end subroutine set_lsf
310
311 real(dp) function get_lsf(rr)
312 real(dp), intent(in) :: rr(NDIM)
313 real(dp) :: qq(NDIM)
314#if NDIM > 1
315 real(dp) :: dist1, dist2
316#endif
317
318 qq = (rr - solution_r0)/solution_radius
319
320#if NDIM == 3
321
322 qq(1) = norm2(qq([1, 2]))
323 qq(2) = qq(3)
324 qq(3) = 0.0_dp
325#endif
326
327 select case (shape)
328 case (1)
329 get_lsf = norm2(qq) - 1.0_dp
330#if NDIM > 1
331 case (2)
332
333 get_lsf = (qq(1)**2 + qq(2)**2 - 1)**3 - &
334 qq(1)**2 * qq(2)**3
335 case (3)
336
337 get_lsf = abs(qq(1))**(2.0_dp/3) / 0.8_dp + &
338 abs(qq(2))**(2.0_dp/3) / 1.5_dp - 0.8_dp
339 case (4)
340
341
342
343 get_lsf = sharpness_t * abs(qq(1)) + abs(qq(2)) - 1.5_dp
344 case (5)
345
346 dist1 = gm_dist_line(rr, [solution_r0(1) - solution_r0(2)/sharpness_t, 0.0_dp], &
347 solution_r0, 2)
348 dist2 = gm_dist_line(rr, [solution_r0(1) + solution_r0(2)/sharpness_t, 0.0_dp], &
349 solution_r0, 2)
350
351
352 qq = rr - solution_r0
353 get_lsf = sharpness_t * abs(qq(1)) + qq(2)
354
355
356 get_lsf = sign(min(dist1, dist2), get_lsf)
357 case (6)
358
359 get_lsf = qq(1)**2 + (qq(2) - abs(qq(1))**(2/3.0_dp))**2 - 1
360 case (7)
361
362 get_lsf = sqrt(sharpness_t * qq(1)**2 + qq(2)**2) - 1.0_dp
363#endif
364 case default
365 error stop "Invalid shape"
366 end select
367
368 end function get_lsf
369
370 subroutine set_error(box)
371 type(box_t), intent(inout) :: box
372 integer :: IJK, nc
373 real(dp) :: rr(NDIM)
374
375 nc = box%n_cell
376 do kji_do(0,nc+1)
377 rr = af_r_cc(box, [ijk])
378 box%cc(ijk, i_error) = box%cc(ijk, i_phi) - solution(rr)
379 end do; close_do
380 end subroutine set_error
381
382 subroutine set_lsf_mask(box)
383 type(box_t), intent(inout) :: box
384 integer :: nc, ix
385
386 nc = box%n_cell
387 box%cc(dtimes(:), i_lsf_mask) = 0
388
389 ix = af_stencil_index(box, mg%operator_key)
390
391 if (ix == af_stencil_none) error stop "No stencil stored"
392
393 if (allocated(box%stencils(ix)%f)) then
394 where (abs(box%stencils(ix)%f) > 0)
395 box%cc(dtimes(1:nc), i_lsf_mask) = 1
396 end where
397 end if
398
399
400
401
402
403
404
405
406
407 end subroutine set_lsf_mask
408
409
410 subroutine bc_solution(box, nb, iv, coords, bc_val, bc_type)
411 type(box_t), intent(in) :: box
412 integer, intent(in) :: nb
413 integer, intent(in) :: iv
414 real(dp), intent(in) :: coords(NDIM, box%n_cell**(NDIM-1))
415 real(dp), intent(out) :: bc_val(box%n_cell**(NDIM-1))
416 integer, intent(out) :: bc_type
417 integer :: n
418
419 bc_type = af_bc_dirichlet
420
421 do n = 1, box%n_cell**(ndim-1)
422 bc_val(n) = solution(coords(:, n))
423 end do
424 end subroutine bc_solution
425
426#if NDIM > 1
427
428
429 subroutine gm_dist_vec_line(r, r0, r1, n_dim, dist_vec, frac)
430 integer, intent(in) :: n_dim
431 real(dp), intent(in) :: r(n_dim), r0(n_dim), r1(n_dim)
432 real(dp), intent(out) :: dist_vec(n_dim)
433 real(dp), intent(out) :: frac
434 real(dp) :: line_len2
435
436 line_len2 = sum((r1 - r0)**2)
437 frac = sum((r - r0) * (r1 - r0))
438
439 if (frac <= 0.0_dp) then
440 frac = 0.0_dp
441 dist_vec = r - r0
442 else if (frac >= line_len2) then
443 frac = 1.0_dp
444 dist_vec = r - r1
445 else
446 dist_vec = r - (r0 + frac/line_len2 * (r1 - r0))
447 frac = sqrt(frac / line_len2)
448 end if
449 end subroutine gm_dist_vec_line
450
451 function gm_dist_line(r, r0, r1, n_dim) result(dist)
452 integer, intent(in) :: n_dim
453 real(dp), intent(in) :: r(n_dim), r0(n_dim), r1(n_dim)
454 real(dp) :: dist, dist_vec(n_dim), frac
455 call gm_dist_vec_line(r, r0, r1, n_dim, dist_vec, frac)
456 dist = norm2(dist_vec)
457 end function gm_dist_line
458#endif
459
460end program poisson_lsf_test
Interface to get variables from the configuration.
Module which contains all Afivo modules, so that a user does not have to include them separately.
Module that allows working with a configuration file.