Octopus
maxwell_boundary_op.F90
Go to the documentation of this file.
1!! Copyright (C) 2019 R. Jestaedt, H. Appel, F. Bonafe, M. Oliveira, N. Tancogne-Dejean
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17
18#include "global.h"
19
21 use accel_oct_m
25 use debug_oct_m
28 use global_oct_m
29 use grid_oct_m
30 use index_oct_m
31 use io_oct_m
33 use math_oct_m
35 use math_oct_m
38 use mesh_oct_m
40 use mpi_oct_m
42 use parser_oct_m
45 use string_oct_m
46 use types_oct_m
47 use unit_oct_m
50 use space_oct_m
52
53 implicit none
54
55 private
56 public :: &
59 bc_mxll_t, &
64
65 type pml_t
66 real(real64) :: width
67 integer :: points_number
68 integer, allocatable :: points_map(:)
69 integer, allocatable :: points_map_inv(:)
70 real(real64) :: power
71 real(real64) :: refl_error
72 real(real64), allocatable :: kappa(:,:)
73 real(real64), allocatable :: sigma_e(:,:)
74 real(real64), allocatable :: sigma_m(:,:)
75 real(real64), allocatable :: a(:,:)
76 real(real64), allocatable :: b(:,:)
77 real(real64), allocatable :: c(:,:)
78 real(real64), allocatable :: mask(:)
79 complex(real64), allocatable :: aux_ep(:,:,:)
80 complex(real64), allocatable :: aux_mu(:,:,:)
81 complex(real64), allocatable :: conv_plus(:,:,:)
82 complex(real64), allocatable :: conv_minus(:,:,:)
83 complex(real64), allocatable :: conv_plus_old(:,:,:)
84 complex(real64), allocatable :: conv_minus_old(:,:,:)
85 logical :: parameters_initialized = .false.
86 ! GPU buffers
87 type(accel_mem_t) :: buff_a, buff_b, buff_c, buff_conv_plus, buff_conv_minus, buff_map, buff_conv_plus_old
88 end type pml_t
89
90 type bc_mxll_t
91 integer :: bc_type(3)
92 integer :: bc_ab_type(3)
93 real(real64) :: bc_bounds(2, 3)
94 logical :: ab_user_def
95 real(real64), allocatable :: ab_ufn(:)
96
97 real(real64) :: ab_width
98 real(real64) :: mask_width
99 integer :: mask_points_number(3)
100 integer, allocatable :: mask_points_map(:,:)
101 real(real64), allocatable :: mask(:,:)
102
103 type(pml_t) :: pml
104 type(single_medium_box_t) :: medium(3)
105
106 integer :: constant_points_number
107 integer, allocatable :: constant_points_map(:)
108 complex(real64), allocatable :: constant_rs_state(:,:)
109 type(accel_mem_t) :: buff_constant_points_map
110
111 integer :: mirror_points_number(3)
112 integer, allocatable :: mirror_points_map(:,:)
113
114 logical :: do_plane_waves = .false. !! look here afterwards
115 type(external_waves_t) :: plane_wave
116 logical :: plane_waves_dims(1:3) = .false.
117
118 real(real64) :: zero_width
119 integer :: zero_points_number(3)
120 integer, allocatable :: zero_points_map(:,:)
121 real(real64), allocatable :: zero(:,:)
122 end type bc_mxll_t
123
124 integer, public, parameter :: &
125 MXLL_BC_ZERO = 0, &
126 mxll_bc_constant = 1, &
127 mxll_bc_mirror_pec = 2, &
128 mxll_bc_mirror_pmc = 3, &
131
132 integer, parameter :: &
133 MXLL_PLANE_WAVES_NEGATIVE_SIDE = -1, &
135
136 integer, public, parameter :: &
137 MXLL_AB_NOT_ABSORBING = 0, &
138 mxll_ab_mask = 1, &
139 mxll_ab_cpml = 2, &
141
142contains
143
144 ! ---------------------------------------------------------
145 subroutine bc_mxll_init(bc, namespace, space, gr, st)
146 type(bc_mxll_t), intent(inout) :: bc
147 type(namespace_t), intent(in) :: namespace
148 class(space_t), intent(in) :: space
149 type(grid_t), intent(in) :: gr
150 type(states_mxll_t), intent(inout) :: st
151
152 integer :: idim, nlines, icol, ncols, ab_shape_dim
153 real(real64) :: bounds(2, space%dim), ab_bounds(2, space%dim)
154 type(block_t) :: blk
155 character(len=1024) :: string
156 character(len=50) :: ab_type_str
157 character(len=1), dimension(3), parameter :: dims = ["x", "y", "z"]
158 logical :: plane_waves_check, ab_mask_check, ab_pml_check, plane_waves_set
159 logical :: constant_check, zero_check
160 real(real64) :: ep_factor, mu_factor, sigma_e_factor, sigma_m_factor
162 push_sub(bc_mxll_init)
164 call profiling_in('BC_MXLL_INIT')
166 plane_waves_set = .false.
167 bc%bc_type(:) = mxll_bc_zero ! default boundary condition is zero
169 !%Variable MaxwellBoundaryConditions
170 !%Type block
171 !%Section Maxwell
172 !%Description
173 !% Defines boundary conditions for the electromagnetic field propagation.
174 !%
175 !% Example:
176 !%
177 !% <tt>%MaxwellBoundaryConditions
178 !% <br>&nbsp;&nbsp; zero | mirror_pec | constant
179 !% <br>%</tt>
180 !%
181 !%
182 !%Option zero 0
183 !% Boundaries are not touched by the Maxwell code - they are set by the normal boundary conditions.
184 !% This is zero boundary conditions for non-periodic systems and periodic boundary conditions if
185 !% PeriodicDimensions is set.
186 !%Option constant 1
187 !% Boundaries are set to a constant.
188 !%Option mirror_pec 2
189 !% Perfect electric conductor.
190 !%Option mirror_pmc 3
191 !% Perfect magnetic conductor.
192 !%Option plane_waves 4
193 !% Boundaries feed in plane waves.
194 !%Option medium 6
195 !% Boundaries as linear medium (not yet implemented).
196 !%End
197 if (parse_block(namespace, 'MaxwellBoundaryConditions', blk) == 0) then
199 call messages_print_with_emphasis(msg='Maxwell Boundary Conditions:', namespace=namespace)
200 ! find out how many lines (i.e. states) the block has
201 nlines = parse_block_n(blk)
202 if (nlines /= 1) then
203 call messages_input_error(namespace, 'MaxwellBoundaryConditions', 'should consist of one line')
204 end if
205 ncols = parse_block_cols(blk, 0)
206 if (ncols /= 3) then
207 call messages_input_error(namespace, 'MaxwellBoundaryConditions', 'should consist of three columns')
208 end if
209 do icol = 1, ncols
210 call parse_block_integer(blk, 0, icol-1, bc%bc_type(icol))
211 end do
212 call parse_block_end(blk)
213 end if
215 do icol = 1, 3
216 select case (bc%bc_type(icol))
217 case (mxll_bc_zero)
218 string = 'Zero'
220 string = 'Constant'
221 case (mxll_bc_mirror_pec)
222 string = 'PEC Mirror'
223 case (mxll_bc_mirror_pmc)
224 string = 'PMC Mirror'
226 string = 'Plane waves'
228 string = 'Medium boundary'
229 case default
230 write(message(1),'(a)') 'Unknown Maxwell boundary condition'
231 call messages_fatal(1, namespace=namespace)
232 end select
233 write(message(1),'(a,I1,a,a)') 'Maxwell boundary condition in direction ', icol, ': ', trim(string)
234 call messages_info(1, namespace=namespace)
235 if (plane_waves_set .and. .not. (parse_is_defined(namespace, 'MaxwellIncidentWaves'))) then
236 write(message(1),'(a)') 'Input: Maxwell boundary condition option is set to "plane_waves".'
237 write(message(2),'(a)') 'Input: User defined Maxwell plane waves have to be defined!'
238 call messages_fatal(2, namespace=namespace)
239 end if
240 end do
241
242 plane_waves_check = .false.
243 ab_mask_check = .false.
244 ab_pml_check = .false.
245 constant_check = .false.
246 zero_check = .false.
247
248 bc%ab_user_def = .false.
249 bc%bc_ab_type(:) = mxll_ab_not_absorbing ! default option
250
251 !%Variable MaxwellAbsorbingBoundaries
252 !%Type block
253 !%Section Maxwell
254 !%Description
255 !% Type of absorbing boundaries used for Maxwell propagation in each direction.
256 !%
257 !% Example:
258 !%
259 !% <tt>%MaxwellAbsorbingBoundaries
260 !% <br>&nbsp;&nbsp; cpml | cpml | cpml
261 !% <br>%</tt>
262 !%
263 !%Option not_absorbing 0
264 !% No absorbing boundaries.
265 !%Option mask 1
266 !% A mask equal to the wavefunctions mask is applied to the Maxwell states at the boundaries
267 !%Option cpml 2
268 !% Perfectly matched layer absorbing boundary
269 !%Option mask_zero 7
270 !% Absorbing boundary region is set to zero
271 !%End
272
273 if (parse_block(namespace, 'MaxwellAbsorbingBoundaries', blk) == 0) then
274 ! find out how many lines (i.e. states) the block has
275 nlines = parse_block_n(blk)
276 if (nlines /= 1) then
277 message(1) = 'MaxwellAbsorbingBounaries has to consist of one line!'
278 call messages_fatal(1, namespace=namespace)
279 end if
280
281 ncols = parse_block_cols(blk, 0)
282 if (ncols /= 3) then
283 message(1) = 'MaxwellAbsorbingBoundaries has to consist of three columns!'
284 call messages_fatal(1, namespace=namespace)
285 end if
286
287 do icol = 1, ncols
288 call parse_block_integer(blk, 0, icol-1, bc%bc_ab_type(icol))
289 end do
290
291 call parse_block_end(blk)
292 end if
293
294 do idim = 1, 3
295 if (bc%bc_ab_type(idim) == mxll_ab_mask) ab_mask_check = .true.
296 if (bc%bc_ab_type(idim) == mxll_ab_cpml) ab_pml_check = .true.
297 if (bc%bc_type(idim) == mxll_bc_constant) constant_check = .true.
298 if (bc%bc_type(idim) == mxll_bc_zero) zero_check = .true.
299 end do
300
301 if (ab_mask_check .or. ab_pml_check) then
302
303 call messages_print_with_emphasis(msg='Maxwell Absorbing Boundaries', namespace=namespace)
304 write(message(1),'(a)') "Please keep in mind that"
305 write(message(2),'(a)') "with larger ABWidth, comes great responsibility."
306 write(message(3),'(a)') "AbsorbingBoundaries occupy space in the simulation box,"
307 write(message(4),'(a)') "hence choose your Lsize wisely."
308 call messages_info(4, namespace=namespace)
309
310 end if
311
312 bounds(1, :) = (gr%idx%nr(2, :) - gr%idx%enlarge(:))*gr%spacing(:)
313 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
314 do idim = 1, st%dim
315 select case (bc%bc_type(idim))
316
317 case (mxll_bc_zero, mxll_bc_mirror_pec, mxll_bc_mirror_pmc)
318
319 bounds(1, idim) = (gr%idx%nr(2, idim) - gr%idx%enlarge(idim))*gr%spacing(idim)
320 bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
321
322 case (mxll_bc_constant)
323
324 bounds(1, idim) = (gr%idx%nr(2, idim) - 2*gr%idx%enlarge(idim))*gr%spacing(idim)
325 bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
326
328
329 bounds(1, idim) = (gr%idx%nr(2, idim) - 2*gr%idx%enlarge(idim))*gr%spacing(idim)
330 bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
331 plane_waves_check = .true.
332 bc%do_plane_waves = .true.
333
334 case (mxll_bc_medium)
335 call bc_mxll_medium_init(gr, namespace, bounds, idim, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
336 call maxwell_medium_points_mapping(bc, gr, bounds)
337 call bc_mxll_generate_medium(bc, space, gr, bounds, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
338
339 end select
340
341 select type (box => gr%box)
342 type is (box_sphere_t)
343 ab_shape_dim = 1
344 if (space%is_periodic()) then
345 message(1) = "Sphere box shape can only work for non-periodic systems"
346 call messages_fatal(1, namespace=namespace)
347 end if
348 type is (box_parallelepiped_t)
349 if (.not. box%axes%orthogonal) then
350 message(1) = "Maxwell propagation does not work for non-orthogonal grids."
351 call messages_fatal(1, namespace=namespace)
352 end if
353
354 ab_shape_dim = space%dim
355 ab_bounds(1, idim) = bounds(1, idim)
356 ab_bounds(2, idim) = bounds(1, idim)
357 class default
358 message(1) = "Box shape for Maxwell propagation not supported yet"
359 call messages_fatal(1, namespace=namespace)
360 end select
361
362 if (bc%bc_ab_type(idim) /= mxll_ab_not_absorbing) then
363
364 call messages_print_var_option("MaxwellAbsorbingBoundaries", bc%bc_ab_type(idim), namespace=namespace)
365
366 select case (bc%bc_ab_type(idim))
367 case (mxll_ab_mask_zero)
368 call bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
369 bc%zero_width = bc%ab_width
370
371 case (mxll_ab_mask)
372 call bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
373 bc%mask_width = bc%ab_width
374
375 case (mxll_ab_cpml)
376 call bc_mxll_pml_init(bc, gr, namespace, ab_bounds, idim)
377
378 case default
379 message(1) = "Absorbing boundary type not implemented for Maxwell propagation"
380 call messages_fatal(1, namespace=namespace)
381 end select
382
383 end if
384
385 select case (bc%bc_ab_type(idim))
387 bounds(1, idim) = ab_bounds(1, idim)
388 bounds(2, idim) = bounds(2, idim)
389 bc%bc_bounds(:, idim) = bounds(:, idim)
390 case default
391 bc%bc_bounds(:, idim) = bounds(:, idim)
392 end select
393
394 select type (box => gr%box)
395 type is (box_parallelepiped_t)
396 select case (bc%bc_ab_type(idim))
397 case (mxll_ab_cpml)
398 ab_type_str = "PML"
399 case (mxll_ab_mask)
400 ab_type_str = "Mask"
401 case (mxll_ab_mask_zero)
402 ab_type_str = "Zero"
403 case default
404 ab_type_str = ""
405 end select
406
407 if (bc%bc_ab_type(idim) == mxll_ab_cpml .or. bc%bc_ab_type(idim) == mxll_ab_mask .or. &
408 bc%bc_ab_type(idim) == mxll_ab_mask_zero) then
409 string = trim(ab_type_str)//" Absorbing Boundary"
410 write(string,'(2a, f10.3,3a)') trim(string), " in "//dims(idim)//' direction spans from: ', &
411 units_from_atomic(units_inp%length, ab_bounds(1, idim) ), ' [', &
412 trim(units_abbrev(units_inp%length)), ']'
413 write(message(1),'(a)') trim(string)
414
415 string = "to "
416 write(string,'(a,f10.3,3a)') trim(string),&
417 units_from_atomic(units_inp%length, ab_bounds(2, idim) ), ' [', &
418 trim(units_abbrev(units_inp%length)), ']'
419
420 write(message(2),'(a)') trim(string)
421 call messages_info(2, namespace=namespace)
422 end if
423
424 class default
425
426 write(message(1),'(a,es10.3,3a)') &
427 " Lower bound = ", units_from_atomic(units_inp%length, ab_bounds(1, idim) ),&
428 ' [', trim(units_abbrev(units_inp%length)), ']'
429 write(message(2),'(a,es10.3,3a)') &
430 " Upper bound = ", units_from_atomic(units_inp%length, ab_bounds(2, idim) ),&
431 ' [', trim(units_abbrev(units_inp%length)), ']'
432 call messages_info(2, namespace=namespace)
433
434 end select
435
436 end do
437
438 ! initialization of surfaces
439 call maxwell_surfaces_init(gr, st, bounds)
440
441 ! mapping of mask boundary points
442 if (ab_mask_check) then
443 call maxwell_mask_points_mapping(bc, gr, ab_bounds)
444 end if
445
446 ! mapping of pml boundary points
447 if (ab_pml_check) then
448 call maxwell_pml_points_mapping(bc, gr, ab_bounds)
449 end if
450
451 ! mapping of constant boundary points
452 if (constant_check) then
453 bounds(1, :) = gr%idx%nr(2, :) * gr%spacing(:)
454 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
455 do idim = 1, st%dim
456 if (bc%bc_type(idim) == mxll_bc_constant) then
457 bounds(1, idim) = (gr%idx%nr(2, idim) - 2*gr%idx%enlarge(idim))*gr%spacing(idim)
458 bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
459 end if
460 end do
461 call maxwell_constant_points_mapping(bc, gr, bounds)
462 end if
463
464 ! mapping of plane waves boundary points
465 if (plane_waves_check) then
466 bounds(1, :) = gr%idx%nr(2, :) * gr%spacing(:)
467 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
468 do idim = 1, st%dim
469 if (bc%bc_type(idim) == mxll_bc_plane_waves) then
470 bounds(1, idim) = (gr%idx%nr(2, idim) - 2*gr%idx%enlarge(idim))*gr%spacing(idim)
471 bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
472 end if
473 end do
474 call maxwell_plane_waves_points_mapping(bc, gr, bounds)
475 call external_waves_init(bc%plane_wave, namespace)
476 end if
477
478 ! mapping of zero points
479 if (zero_check) then
480 bounds(1, :) = gr%idx%nr(2, :) * gr%spacing(:)
481 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
482 do idim = 1, st%dim
483 if (bc%bc_type(idim) == mxll_bc_zero) then
484 bounds(1, idim) = (gr%idx%nr(2, idim) - gr%idx%enlarge(idim))*gr%spacing(idim)
485 bounds(2, idim) = (gr%idx%nr(2, idim)) * gr%spacing(idim)
486 end if
487 end do
488 call maxwell_zero_points_mapping(bc, gr, bounds)
489 end if
490
491 if (ab_mask_check) then
492 call bc_mxll_generate_mask(bc, gr, ab_bounds)
493 end if
494
495 if (ab_pml_check) then
496 call bc_mxll_generate_pml(bc, space)
497 end if
498
499 !call bc_generate_zero(bc, gr, ab_bounds)
500
501 if (debug%info) call bc_mxll_write_info(bc, gr, namespace, space)
502
503 if (ab_mask_check .or. ab_pml_check) then
504 call messages_print_with_emphasis(namespace=namespace)
505 end if
506
507 call profiling_out('BC_MXLL_INIT')
508
509 pop_sub(bc_mxll_init)
510 end subroutine bc_mxll_init
511
512 ! ---------------------------------------------------------
513 subroutine bc_mxll_end(bc)
514 type(bc_mxll_t), intent(inout) :: bc
515
516 integer :: idim
517
518 push_sub(bc_mxll_end)
519
520 safe_deallocate_a(bc%ab_ufn)
521
522 safe_deallocate_a(bc%mask_points_map)
523 safe_deallocate_a(bc%mask)
524
525 call pml_end(bc%pml)
526 do idim = 1, 3
527 call single_medium_box_end(bc%medium(idim))
528 end do
529
530 safe_deallocate_a(bc%constant_points_map)
531 safe_deallocate_a(bc%constant_rs_state)
532 if (accel_is_enabled()) then
533 call accel_release_buffer(bc%buff_constant_points_map)
534 end if
535
536 safe_deallocate_a(bc%mirror_points_map)
537
538 call external_waves_end(bc%plane_wave)
539
540 safe_deallocate_a(bc%zero_points_map)
541 safe_deallocate_a(bc%zero)
542
543 pop_sub(bc_mxll_end)
544 end subroutine bc_mxll_end
545
546 ! ---------------------------------------------------------
547 subroutine pml_end(pml)
548 type(pml_t), intent(inout) :: pml
549
550 push_sub(pml_end)
551
552 safe_deallocate_a(pml%points_map)
553 safe_deallocate_a(pml%points_map_inv)
554 safe_deallocate_a(pml%kappa)
555 safe_deallocate_a(pml%sigma_e)
556 safe_deallocate_a(pml%sigma_m)
557 safe_deallocate_a(pml%a)
558 safe_deallocate_a(pml%b)
559 safe_deallocate_a(pml%c)
560 safe_deallocate_a(pml%mask)
561 safe_deallocate_a(pml%aux_ep)
562 safe_deallocate_a(pml%aux_mu)
563 safe_deallocate_a(pml%conv_plus)
564 safe_deallocate_a(pml%conv_minus)
565 safe_deallocate_a(pml%conv_plus_old)
566 safe_deallocate_a(pml%conv_minus_old)
567 if (accel_is_enabled()) then
568 call accel_release_buffer(pml%buff_map)
569 call accel_release_buffer(pml%buff_a)
570 call accel_release_buffer(pml%buff_b)
571 call accel_release_buffer(pml%buff_c)
572 call accel_release_buffer(pml%buff_conv_plus)
573 call accel_release_buffer(pml%buff_conv_minus)
574 call accel_release_buffer(pml%buff_conv_plus_old)
575 end if
576
577
578 pop_sub(pml_end)
579 end subroutine pml_end
580
581
582
583 ! ---------------------------------------------------------
584 subroutine bc_mxll_medium_init(gr, namespace, bounds, idim, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
585 type(grid_t), intent(in) :: gr
586 type(namespace_t), intent(in) :: namespace
587 real(real64), intent(inout) :: bounds(:,:)
588 integer, intent(in) :: idim
589 real(real64), intent(out) :: ep_factor
590 real(real64), intent(out) :: mu_factor
591 real(real64), intent(out) :: sigma_e_factor
592 real(real64), intent(out) :: sigma_m_factor
593
594 real(real64) :: width
595
596 push_sub(bc_mxll_medium_init)
597
598 call profiling_in('BC_MXLL_MEDIUM_INIT')
599
600 !%Variable MediumWidth
601 !%Type float
602 !%Default 0.0 a.u.
603 !%Section Maxwell::Boundaries
604 !%Description
605 !% Width of the boundary region with medium
606 !%End
607 call parse_variable(namespace, 'MediumWidth', m_zero, width, units_inp%length)
608 bounds(1,idim) = ( gr%idx%nr(2, idim) - gr%idx%enlarge(idim) ) * gr%spacing(idim)
609 bounds(1,idim) = bounds(1,idim) - width
610 bounds(2,idim) = ( gr%idx%nr(2, idim) ) * gr%spacing(idim)
611
612 !%Variable MediumEpsilonFactor
613 !%Type float
614 !%Default 1.0.
615 !%Section Maxwell::Boundaries
616 !%Description
617 !% Linear medium electric susceptibility.
618 !%End
619 call parse_variable(namespace, 'MediumpsilonFactor', m_one, ep_factor, unit_one)
620
621 !%Variable MediumMuFactor
622 !%Type float
623 !%Default 1.0
624 !%Section Maxwell::Boundaries
625 !%Description
626 !% Linear medium magnetic susceptibility.
627 !%End
628 call parse_variable(namespace, 'MediumMuFactor', m_one, mu_factor, unit_one)
629
630 !%Variable MediumElectricSigma
631 !%Type float
632 !%Default 0.
633 !%Section Maxwell::Boundaries
634 !%Description
635 !% Electric conductivity of the linear medium.
636 !%End
637
638 call parse_variable(namespace, 'MediumElectricSigma', m_zero, sigma_e_factor, unit_one)
639 !%Variable MediumMagneticSigma
640 !%Type float
641 !%Default 0.
642 !%Section Maxwell::Boundaries
643 !%Description
644 !% Magnetic conductivity of the linear medium.
645 !%End
646 call parse_variable(namespace, 'MediumMagneticSigma', m_zero, sigma_m_factor, unit_one)
647
648 call profiling_out('BC_MXLL_MEDIUM_INIT')
649
650 pop_sub(bc_mxll_medium_init)
651 end subroutine bc_mxll_medium_init
652
653 ! ---------------------------------------------------------
654 subroutine bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
655 type(bc_mxll_t), intent(inout) :: bc
656 type(grid_t), intent(in) :: gr
657 type(namespace_t), intent(in) :: namespace
658 real(real64), intent(inout) :: ab_bounds(:,:)
659 integer, intent(in) :: idim
660
661 real(real64) :: default_width
662
663 push_sub(bc_mxll_ab_bounds_init)
664
665 !%Variable MaxwellABWidth
666 !%Type float
667 !%Section Maxwell::Boundaries
668 !%Description
669 !% Width of the region used to apply the absorbing boundaries. The default value is twice
670 !% the derivative order.
671 !%End
672
673 default_width = m_two * gr%der%order * maxval(gr%spacing(1:3))
674 call parse_variable(namespace, 'MaxwellABWidth', default_width, bc%ab_width, units_inp%length)
675
676 if (bc%ab_width < default_width) then
677 bc%ab_width = default_width
678 write(message(1),'(a)') 'Absorbing boundary width has to be larger or equal than twice the derivatives order times spacing!'
679 write(message(2),'(a,es10.3)') 'Absorbing boundary width is set to: ', default_width
680 call messages_info(2, namespace=namespace)
681 end if
682
683 ab_bounds(1, idim) = ab_bounds(2, idim) - bc%ab_width
684
686 end subroutine bc_mxll_ab_bounds_init
687
688 ! ---------------------------------------------------------
689 subroutine bc_mxll_pml_init(bc, gr, namespace, ab_bounds, idim)
690 type(bc_mxll_t), intent(inout) :: bc
691 type(grid_t), intent(in) :: gr
692 type(namespace_t), intent(in) :: namespace
693 real(real64), intent(inout) :: ab_bounds(:,:)
694 integer, intent(in) :: idim
695
696 push_sub(bc_mxll_pml_init)
697
698 call bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
699 bc%pml%width = bc%ab_width
700
701 !%Variable MaxwellABPMLPower
702 !%Type float
703 !%Default 3.5
704 !%Section Maxwell::Boundaries
705 !%Description
706 !% Exponential of the polynomial profile for the non-physical conductivity of the PML.
707 !% Should be between 2 and 4
708 !%End
709 call parse_variable(namespace, 'MaxwellABPMLPower', 3.5_real64, bc%pml%power, unit_one)
710
711 !%Variable MaxwellABPMLReflectionError
712 !%Type float
713 !%Default 1.0e-16
714 !%Section Maxwell::Boundaries
715 !%Description
716 !% Tolerated reflection error for the PML
717 !%End
718 call parse_variable(namespace, 'MaxwellABPMLReflectionError', 1.0e-16_real64, bc%pml%refl_error, unit_one)
719
720 pop_sub(bc_mxll_pml_init)
721 end subroutine bc_mxll_pml_init
722
723 ! ---------------------------------------------------------
724 subroutine bc_mxll_write_info(bc, mesh, namespace, space)
725 type(bc_mxll_t), intent(in) :: bc
726 class(mesh_t), intent(in) :: mesh
727 type(namespace_t), intent(in) :: namespace
728 class(space_t), intent(in) :: space
729
730 integer :: err, idim, idim2
731 real(real64), allocatable :: tmp(:)
732 logical :: mask_check, pml_check, medium_check
733 character(1) :: dim_label(3)
734
735 push_sub(bc_mxll_write_info)
736
737 call profiling_in('BC_MXLL_WRITE_INFO')
738
739 mask_check = .false.
740 pml_check = .false.
741 medium_check = .false.
742
743 dim_label = (/'x', 'y', 'z'/)
744
745 do idim = 1, 3
746 if (bc%bc_ab_type(idim) == mxll_ab_mask) then
747 mask_check = .true.
748 end if
749 if (bc%bc_ab_type(idim) == mxll_ab_cpml) then
750 pml_check = .true.
751 end if
752 if (bc%bc_ab_type(idim) == mxll_bc_medium) then
753 medium_check = .true.
754 end if
755 end do
756
757 if (mask_check) then
758 safe_allocate(tmp(1:mesh%np))
759 do idim = 1, 3
760 tmp(:) = m_zero
761 call get_mask_io_function(bc%mask, bc, tmp, idim)
762 call write_files("maxwell_mask", tmp)
763 end do
764 safe_deallocate_a(tmp)
765 else if (pml_check) then
766 safe_allocate(tmp(1:mesh%np))
767 do idim = 1, 3
768 ! sigma for electric field dim = idim
769 tmp(:) = m_one
770 call get_pml_io_function(bc%pml%sigma_e(:, idim), bc, tmp)
771 call write_files("maxwell_sigma_e-"//dim_label(idim), tmp)
772
773 ! sigma for magnetic dim = idim
774 tmp(:) = m_zero
775 call get_pml_io_function(bc%pml%sigma_m(:, 1), bc, tmp)
776 call write_files("maxwell_sigma_m-"//dim_label(idim), tmp)
777
778 ! pml_a for electric field dim = idim
779 tmp(:) = m_zero
780 call get_pml_io_function(real(bc%pml%a(:, idim), real64), bc, tmp)
781 call write_files("maxwell_sigma_pml_a_e-"//dim_label(idim), tmp)
782 end do
783 safe_deallocate_a(tmp)
784 end if
785
786 if (medium_check) then
787 safe_allocate(tmp(1:mesh%np))
788 do idim = 1, 3
789 ! medium epsilon
790 tmp(:) = p_ep
791 call get_medium_io_function(bc%medium(idim)%ep, bc, tmp, idim)
792 call write_files("maxwell_ep"//dim_label(idim), tmp)
793 ! medium mu
794 tmp(:) = p_mu
795 call get_medium_io_function(bc%medium(idim)%mu, bc, tmp, idim)
796 call write_files("maxwell_mu"//dim_label(idim), tmp)
797 ! medium epsilon
798 tmp(:) = p_c
799 call get_medium_io_function(bc%medium(idim)%c, bc, tmp, idim)
800 call write_files("maxwell_c"//dim_label(idim), tmp)
801 do idim2 = 1, 3
802 ! medium epsilon aux field dim = idim
803 tmp(:) = m_zero
804 call get_medium_io_function(bc%medium(idim)%aux_ep(:, idim2), bc, tmp, idim)
805 call write_files("maxwell_aux_ep-"//dim_label(idim)//"-"//dim_label(idim2), tmp)
806
807 ! medium mu aux field dim = idim
808 tmp(:) = m_zero
809 call get_medium_io_function(bc%medium(idim)%aux_mu(:, idim2), bc, tmp, idim)
810 call write_files("maxwell_aux_mu-"//dim_label(idim)//"-"//dim_label(idim2), tmp)
811 end do
812 end do
813
814 safe_deallocate_a(tmp)
815 end if
816
817 call profiling_out('BC_MXLL_WRITE_INFO')
818
820 contains
821
822 subroutine get_pml_io_function(pml_func, bc, io_func)
823 real(real64), intent(in) :: pml_func(:)
824 type(bc_mxll_t), intent(in) :: bc
825 real(real64), intent(inout) :: io_func(:)
826
827 integer :: ip, ip_in
828
829 do ip_in = 1, bc%pml%points_number
830 ip = bc%pml%points_map(ip_in)
831 io_func(ip) = pml_func(ip_in)
832 end do
833
834 end subroutine get_pml_io_function
835
836 subroutine get_mask_io_function(mask_func, bc, io_func, idim)
837 real(real64), intent(in) :: mask_func(:,:)
838 type(bc_mxll_t), intent(in) :: bc
839 real(real64), intent(inout) :: io_func(:)
840 integer, intent(in) :: idim
841
842 integer :: ip, ip_in
843
844 do ip_in = 1, bc%mask_points_number(idim)
845 ip = bc%mask_points_map(ip_in, idim)
846 io_func(ip) = mask_func(ip_in, idim)
847 end do
848
849 end subroutine get_mask_io_function
850
851 subroutine get_medium_io_function(medium_func, bc, io_func, idim)
852 real(real64), intent(in) :: medium_func(:)
853 type(bc_mxll_t), intent(in) :: bc
854 real(real64), intent(inout) :: io_func(:)
855 integer, intent(in) :: idim
856
857 integer :: ip, ip_in
858
859 do ip_in = 1, bc%medium(idim)%points_number
860 ip = bc%medium(idim)%points_map(ip_in)
861 io_func(ip) = medium_func(ip_in)
862 end do
863
864 end subroutine get_medium_io_function
865
866 subroutine write_files(filename, tmp)
867 character(len=*), intent(in) :: filename
868 real(real64), intent(in) :: tmp(:)
869
870 call dio_function_output(io_function_fill_how("VTK"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
871 unit_one, err)
872 call dio_function_output(io_function_fill_how("AxisX"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
873 unit_one, err)
874 call dio_function_output(io_function_fill_how("AxisY"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
875 unit_one, err)
876 call dio_function_output(io_function_fill_how("AxisZ"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
877 unit_one, err)
878 call dio_function_output(io_function_fill_how("PlaneX"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
879 unit_one, err)
880 call dio_function_output(io_function_fill_how("PlaneY"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
881 unit_one, err)
882 call dio_function_output(io_function_fill_how("PlaneZ"), "./td.general", trim(filename), namespace, space, mesh, tmp, &
883 unit_one, err)
884 end subroutine write_files
885
886 end subroutine bc_mxll_write_info
887
888 ! ---------------------------------------------------------
889 subroutine maxwell_mask_points_mapping(bc, mesh, bounds)
890 type(bc_mxll_t), intent(inout) :: bc
891 class(mesh_t), intent(in) :: mesh
892 real(real64), intent(in) :: bounds(:,:)
893
894 integer :: ip, ip_in, ip_in_max, point_info, idim
895
897
898 call profiling_in('MAXWELL_MASK_POINTS_MAPPING')
899 ip_in_max = 1
900 do idim = 1, 3
901 if (bc%bc_ab_type(idim) == mxll_ab_mask) then
902 ! allocate mask points map
903 ip_in = 0
904 do ip = 1, mesh%np
905 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
906 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
907 ip_in = ip_in + 1
908 end if
909 end do
910 if (ip_in > ip_in_max) ip_in_max = ip_in
911 bc%mask_points_number(idim) = ip_in
912 end if
913 end do
914 safe_allocate(bc%mask(1:ip_in_max, 1:idim))
915 safe_allocate(bc%mask_points_map(1:ip_in_max, 1:idim))
916
917 do idim = 1, 3
918 if (bc%bc_ab_type(idim) == mxll_ab_mask) then
919 ! mask points mapping
920 ip_in = 0
921 do ip = 1, mesh%np
922 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
923 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
924 ip_in = ip_in + 1
925 bc%mask_points_map(ip_in, idim) = ip
926 end if
927 end do
928 end if
929 end do
930
931 call profiling_out('MAXWELL_MASK_POINTS_MAPPING')
933 end subroutine maxwell_mask_points_mapping
934
935 ! ---------------------------------------------------------
936 subroutine maxwell_pml_points_mapping(bc, mesh, bounds)
937 type(bc_mxll_t), intent(inout) :: bc
938 class(mesh_t), intent(in) :: mesh
939 real(real64), intent(in) :: bounds(:,:)
940
941 integer :: ip, ip_in, point_info
942
944
945 call profiling_in('MXWL_PML_POINTS_MAPPING')
947 ! allocate pml points map
948 ip_in = 0
949 do ip = 1, mesh%np
950 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
951 if (point_info == 1) then
952 ip_in = ip_in + 1
953 end if
954 end do
955 bc%pml%points_number = ip_in
956 safe_allocate(bc%pml%points_map(1:ip_in))
957 safe_allocate(bc%pml%points_map_inv(1:mesh%np))
958 bc%pml%points_map_inv = 0
959 ip_in = 0
960 do ip = 1, mesh%np
961 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
962 if (point_info == 1) then
963 ip_in = ip_in + 1
964 bc%pml%points_map(ip_in) = ip
965 bc%pml%points_map_inv(ip) = ip_in
966 end if
967 end do
968
969 call profiling_out('MXWL_PML_POINTS_MAPPING')
970
972 end subroutine maxwell_pml_points_mapping
973
974 ! ---------------------------------------------------------
975 subroutine maxwell_constant_points_mapping(bc, mesh, bounds)
976 type(bc_mxll_t), intent(inout) :: bc
977 class(mesh_t), intent(in) :: mesh
978 real(real64), intent(in) :: bounds(:,:)
979
980 integer :: ip, ip_in, point_info
981
983
984 call profiling_in('MAXWELL_CONSTANT_POINTS_MAP')
985
986 ! allocate constant points map
987 ip_in = 0
988 do ip = 1, mesh%np
989 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
990 if (point_info == 1) then
991 ip_in = ip_in + 1
992 end if
993 end do
994 bc%constant_points_number = ip_in
995 safe_allocate(bc%constant_points_map(1:ip_in))
996 safe_allocate(bc%constant_rs_state(1:ip_in, 1:3))
997
998 ! zero constant mapping
999 ip_in = 0
1000 do ip = 1, mesh%np
1001 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1002 if (point_info == 1) then
1003 ip_in = ip_in + 1
1004 bc%constant_points_map(ip_in) = ip
1005 end if
1006 end do
1007
1008 if (accel_is_enabled()) then
1009 call accel_create_buffer(bc%buff_constant_points_map, accel_mem_read_only, type_integer, bc%constant_points_number)
1010 call accel_write_buffer(bc%buff_constant_points_map, bc%constant_points_number, bc%constant_points_map)
1011 end if
1012
1013 call profiling_out('MAXWELL_CONSTANT_POINTS_MAP')
1014
1016 end subroutine maxwell_constant_points_mapping
1017
1018 ! ---------------------------------------------------------
1019 subroutine maxwell_plane_waves_points_mapping(bc, mesh, bounds)
1020 type(bc_mxll_t), intent(inout) :: bc
1021 class(mesh_t), intent(in) :: mesh
1022 real(real64), intent(in) :: bounds(:,:)
1023
1024 integer :: ip, ip_in, point_info
1025
1027
1028 call profiling_in('MXLL_PW_POINTS_MAP')
1029
1030 bc%plane_waves_dims = (bc%bc_type(1:3) == mxll_bc_plane_waves)
1032 ! allocate map
1033 ip_in = 0
1034 do ip = 1, mesh%np
1035 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1036 if (point_info == 1) then
1037 ip_in = ip_in + 1
1038 end if
1039 end do
1040 bc%plane_wave%points_number = ip_in
1041 safe_allocate(bc%plane_wave%points_map(1:ip_in))
1042
1043 ! mapping
1044 ip_in = 0
1045 do ip = 1, mesh%np
1046 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1047 if (point_info == 1) then
1048 ip_in = ip_in + 1
1049 bc%plane_wave%points_map(ip_in) = ip
1050 end if
1051 end do
1052
1053 if (accel_is_enabled()) then
1054 call accel_create_buffer(bc%plane_wave%buff_map, accel_mem_read_only, type_integer, bc%plane_wave%points_number)
1055 call accel_write_buffer(bc%plane_wave%buff_map, bc%plane_wave%points_number, bc%plane_wave%points_map)
1056 end if
1057
1058 call profiling_out('MXLL_PW_POINTS_MAP')
1059
1062
1063 ! ---------------------------------------------------------
1064 subroutine maxwell_zero_points_mapping(bc, mesh, bounds)
1065 type(bc_mxll_t), intent(inout) :: bc
1066 class(mesh_t), intent(in) :: mesh
1067 real(real64), intent(in) :: bounds(:,:)
1068
1069 integer :: ip, ip_in, ip_in_max, point_info, idim
1072
1073 call profiling_in('MXLL_ZERO_POINTS_MAPPING')
1074
1075 ip_in_max = 0
1076 do idim = 1, 3
1077 if (bc%bc_type(idim) == mxll_bc_zero) then
1078 ! allocate zero points map
1079 ip_in = 0
1080 do ip = 1, mesh%np
1081 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1082 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
1083 ip_in = ip_in + 1
1084 end if
1085 end do
1086 if (ip_in > ip_in_max) ip_in_max = ip_in
1087 end if
1088 end do
1089 bc%zero_points_number = ip_in
1090 safe_allocate(bc%zero(1:ip_in_max,3))
1091 safe_allocate(bc%zero_points_map(1:ip_in_max, 1:3))
1092
1093 do idim = 1, 3
1094 if (bc%bc_type(idim) == mxll_bc_zero) then
1095 ! zero points mapping
1096 ip_in = 0
1097 do ip = 1, mesh%np
1098 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1099 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
1100 ip_in = ip_in + 1
1101 bc%zero_points_map(ip_in, idim) = ip
1102 end if
1103 end do
1104 end if
1105 end do
1106
1107 call profiling_out('MXLL_ZERO_POINTS_MAPPING')
1108
1110 end subroutine maxwell_zero_points_mapping
1111
1112 ! ---------------------------------------------------------
1113 subroutine maxwell_medium_points_mapping(bc, mesh, bounds)
1114 type(bc_mxll_t), intent(inout) :: bc
1115 class(mesh_t), intent(in) :: mesh
1116 real(real64), intent(in) :: bounds(:,:)
1117
1118 integer :: ip, ip_in, ip_in_max, ip_bd, ip_bd_max, point_info, boundary_info, idim
1119
1121
1122 call profiling_in('MXLL_MEDIUM_POINTS_MAPPING')
1123
1124 ip_in_max = 0
1125 ip_bd_max = 0
1126 do idim = 1, 3
1127 if (bc%bc_type(idim) == mxll_bc_medium) then
1128 ! allocate pml points map
1129 ip_in = 0
1130 ip_bd = 0
1131 do ip = 1, mesh%np
1132 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1133 call maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
1134 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
1135 ip_in = ip_in + 1
1136 end if
1137 if ((boundary_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
1138 ip_bd = ip_bd + 1
1139 end if
1140 end do
1141 bc%medium(idim)%points_number = ip_in
1142 bc%medium(idim)%bdry_number = ip_bd
1143 end if
1144 end do
1145 do idim = 1, 3
1146 safe_allocate(bc%medium(idim)%points_map(1:ip_in_max))
1147 safe_allocate(bc%medium(idim)%bdry_map(1:ip_bd_max))
1148 end do
1149
1150 ip_in = 0
1151 ip_bd = 0
1152 do idim = 1, 3
1153 if (bc%bc_type(idim) == mxll_bc_medium) then
1154 do ip = 1, mesh%np
1155 call maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1156 call maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
1157 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
1158 ip_in = ip_in + 1
1159 bc%medium(idim)%points_map(ip_in) = ip
1160 end if
1161 if ((boundary_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim))) then
1162 ip_bd = ip_bd + 1
1163 bc%medium(idim)%bdry_map(ip_bd) = ip
1164 end if
1165 end do
1166 end if
1167 end do
1168
1169 call profiling_out('MXLL_MEDIUM_POINTS_MAPPING')
1170
1172 end subroutine maxwell_medium_points_mapping
1173
1174 ! ---------------------------------------------------------
1175 subroutine bc_mxll_generate_pml(bc, space)
1176 type(bc_mxll_t), intent(inout) :: bc
1177 class(space_t), intent(in) :: space
1178
1179
1180 push_sub(bc_mxll_generate_pml)
1181
1182 call profiling_in('BC_MXLL_GENERATE_PML')
1183
1184 safe_allocate(bc%pml%kappa(1:bc%pml%points_number, 1:space%dim))
1185 safe_allocate(bc%pml%sigma_e(1:bc%pml%points_number, 1:space%dim))
1186 safe_allocate(bc%pml%sigma_m(1:bc%pml%points_number, 1:space%dim))
1187 safe_allocate(bc%pml%a(1:bc%pml%points_number, 1:space%dim))
1188 safe_allocate(bc%pml%b(1:bc%pml%points_number, 1:space%dim))
1189 safe_allocate(bc%pml%c(1:bc%pml%points_number, 1:3))
1190 safe_allocate(bc%pml%mask(1:bc%pml%points_number))
1191 safe_allocate(bc%pml%conv_plus(1:bc%pml%points_number, 1:3, 1:3))
1192 safe_allocate(bc%pml%conv_minus(1:bc%pml%points_number, 1:3, 1:3))
1193 safe_allocate(bc%pml%conv_plus_old(1:bc%pml%points_number, 1:3, 1:3))
1194 safe_allocate(bc%pml%conv_minus_old(1:bc%pml%points_number, 1:3, 1:3))
1195 safe_allocate(bc%pml%aux_ep(1:bc%pml%points_number, 1:3, 1:3))
1196 safe_allocate(bc%pml%aux_mu(1:bc%pml%points_number, 1:3, 1:3))
1197
1198 bc%pml%kappa = m_one
1199 bc%pml%sigma_e = m_zero
1200 bc%pml%sigma_m = m_zero
1201 bc%pml%a = m_zero
1202 bc%pml%b = m_zero
1203 bc%pml%c = m_zero
1204 bc%pml%mask = m_one
1205 bc%pml%conv_plus = m_z0
1206 bc%pml%conv_minus = m_z0
1207 bc%pml%conv_plus_old = m_z0
1208 bc%pml%conv_minus_old = m_z0
1209
1210 if (accel_is_enabled()) then
1211 call accel_create_buffer(bc%pml%buff_map, accel_mem_read_only, type_integer, bc%pml%points_number)
1212 call accel_create_buffer(bc%pml%buff_a, accel_mem_read_only, type_float, int(bc%pml%points_number, int64)*space%dim)
1213 call accel_create_buffer(bc%pml%buff_b, accel_mem_read_only, type_float, int(bc%pml%points_number, int64)*space%dim)
1214 call accel_create_buffer(bc%pml%buff_c, accel_mem_read_only, type_float, int(bc%pml%points_number, int64)*space%dim)
1215 call accel_create_buffer(bc%pml%buff_conv_plus, accel_mem_read_write, type_cmplx, &
1216 int(bc%pml%points_number, int64)*space%dim*space%dim)
1217 call accel_create_buffer(bc%pml%buff_conv_minus, accel_mem_read_write, type_cmplx, &
1218 int(bc%pml%points_number, int64)*space%dim*space%dim)
1219 call accel_create_buffer(bc%pml%buff_conv_plus_old, accel_mem_read_write, type_cmplx, &
1220 int(bc%pml%points_number, int64)*space%dim*space%dim)
1221
1222 call accel_write_buffer(bc%pml%buff_a, bc%pml%points_number, space%dim, bc%pml%a)
1223 call accel_write_buffer(bc%pml%buff_b, bc%pml%points_number, space%dim, bc%pml%b)
1224 call accel_write_buffer(bc%pml%buff_c, bc%pml%points_number, space%dim, bc%pml%c)
1225 call accel_write_buffer(bc%pml%buff_conv_plus, bc%pml%points_number, space%dim, space%dim, bc%pml%conv_plus)
1226 call accel_write_buffer(bc%pml%buff_conv_minus, bc%pml%points_number, space%dim, space%dim, bc%pml%conv_minus)
1227 call accel_write_buffer(bc%pml%buff_conv_plus_old, bc%pml%points_number, space%dim, space%dim, bc%pml%conv_plus_old)
1228 end if
1229
1230 call profiling_out('BC_MXLL_GENERATE_PML')
1231
1232 pop_sub(bc_mxll_generate_pml)
1233 end subroutine bc_mxll_generate_pml
1234
1235
1236 ! ---------------------------------------------------------
1237 subroutine bc_mxll_generate_pml_parameters(bc, space, gr, c_factor, dt)
1238 type(bc_mxll_t), intent(inout) :: bc
1239 class(space_t), intent(in) :: space
1240 type(grid_t), intent(in) :: gr
1241 real(real64), intent(in) :: c_factor
1242 real(real64), optional, intent(in) :: dt
1243
1244 integer :: ip, ip_in, idim
1245 real(real64) :: width(3)
1246 real(real64) :: ddv(3), ss_e, ss_m, ss_max, aa_e, aa_m, bb_e, bb_m, gg, hh, kk, ll_e, ll_m
1247 real(real64), allocatable :: tmp(:), tmp_grad(:,:)
1248
1249
1251 safe_allocate(tmp(gr%np_part))
1252 safe_allocate(tmp_grad(gr%np, 1:space%dim))
1253 ! here bounds are ab_bounds, for which ab_bounds(1,:) = bc%bounds(1,:)
1254 ! width is stored in bc%pml%width
1255 ! assuming that the width is the same in the 3 dimensions (only case implemented now), we can change the following line:
1256 ! width(1:3) = bounds(2, 1:3) - bounds(1, 1:3) !! original line
1257 ! by
1258 width(1:3) = (/ bc%pml%width, bc%pml%width, bc%pml%width /)
1259
1260 ! PML variables for all boundary points
1261 do ip_in = 1, bc%pml%points_number
1262 ip = bc%pml%points_map(ip_in)
1263 ddv(1:3) = abs(gr%x(ip, 1:3)) - bc%bc_bounds(1, 1:3)
1264 do idim = 1, space%dim
1265 if (ddv(idim) >= m_zero) then
1266 gg = (ddv(idim)/bc%pml%width)**bc%pml%power
1267 hh = (m_one-ddv(idim)/bc%pml%width)**bc%pml%power
1268 kk = m_one
1269 ss_max = -(bc%pml%power + m_one)*p_c*c_factor*p_ep*log(bc%pml%refl_error)/(m_two * bc%pml%width)
1270 ss_e = gg * ss_max
1271 ss_m = gg * ss_max
1272 ll_e = ss_e*kk
1273 ll_m = ss_m*kk
1274 bb_e = exp(-(ss_e/(p_ep))*dt)
1275 bb_m = exp(-(ss_m/(p_ep))*dt)
1276 aa_e = (bb_e - 1)
1277 aa_m = (bb_m - 1)
1278 if (abs(ll_e) < m_epsilon) aa_e = m_zero
1279 if (abs(ll_m) < m_epsilon) aa_m = m_zero
1280 bc%pml%sigma_e(ip_in, idim) = ss_e
1281 bc%pml%sigma_m(ip_in, idim) = ss_m
1282 bc%pml%a(ip_in, idim) = aa_e
1283 bc%pml%b(ip_in, idim) = bb_e
1284 bc%pml%kappa(ip_in, idim) = kk
1285 bc%pml%mask(ip_in) = bc%pml%mask(ip_in) * (m_one - sin(ddv(idim)*m_pi/(m_two*(width(idim))))**2)
1286 else
1287 bc%pml%kappa(ip_in, idim) = m_one
1288 bc%pml%sigma_e(ip_in, idim) = m_zero
1289 bc%pml%sigma_m(ip_in, idim) = m_zero
1290 bc%pml%a(ip_in, idim) = m_zero
1291 bc%pml%b(ip_in, idim) = m_zero
1292 bc%pml%mask(ip_in) = m_one
1293 end if
1294 end do
1295 end do
1296
1297 ! PML auxiliary epsilon for all boundary points
1298 do idim = 1, space%dim
1299 tmp = p_ep
1300 do ip_in = 1, bc%pml%points_number
1301 ip = bc%pml%points_map(ip_in)
1302 tmp(ip) = p_ep / bc%pml%kappa(ip_in, idim)
1303 end do
1304 call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
1305 do ip_in = 1, bc%pml%points_number
1306 ip = bc%pml%points_map(ip_in)
1307 bc%pml%aux_ep(ip_in, :, idim) = tmp_grad(ip, :)/(m_four*p_ep*bc%pml%kappa(ip_in, idim))
1308 end do
1309 end do
1310
1311 ! PML auxiliary mu
1312 do idim = 1, space%dim
1313 ! P_mu is proportional to 1/P_c**2, so we need to devide by c_factor
1314 tmp = p_mu/c_factor**2
1315 do ip_in = 1, bc%pml%points_number
1316 ip = bc%pml%points_map(ip_in)
1317 tmp(ip) = p_mu/c_factor**2 / bc%pml%kappa(ip_in, idim)
1318 end do
1319 call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
1320 do ip_in = 1, bc%pml%points_number
1321 ip = bc%pml%points_map(ip_in)
1322 bc%pml%aux_mu(ip_in, :, idim) = tmp_grad(ip, :)/(m_four*p_mu/c_factor**2 * bc%pml%kappa(ip_in, idim))
1323 end do
1324 end do
1325
1326 ! PML auxiliary c for all boundary points
1327 do idim = 1, space%dim
1328 do ip_in = 1, bc%pml%points_number
1329 bc%pml%c(ip_in, idim) = p_c*c_factor/bc%pml%kappa(ip_in, idim)
1330 end do
1331 end do
1332 safe_deallocate_a(tmp)
1333 safe_deallocate_a(tmp_grad)
1334
1335 if (accel_is_enabled()) then
1336 call accel_write_buffer(bc%pml%buff_map, bc%pml%points_number, bc%pml%points_map)
1337 call accel_write_buffer(bc%pml%buff_a, bc%pml%points_number, space%dim, bc%pml%a)
1338 call accel_write_buffer(bc%pml%buff_b, bc%pml%points_number, space%dim, bc%pml%b)
1339 call accel_write_buffer(bc%pml%buff_c, bc%pml%points_number, space%dim, bc%pml%c)
1340 call accel_write_buffer(bc%pml%buff_conv_plus, bc%pml%points_number, space%dim, space%dim, bc%pml%conv_plus)
1341 call accel_write_buffer(bc%pml%buff_conv_minus, bc%pml%points_number, space%dim, space%dim, bc%pml%conv_minus)
1342 end if
1343
1344 bc%pml%parameters_initialized = .true.
1345
1347
1348 end subroutine bc_mxll_generate_pml_parameters
1349
1350 ! ---------------------------------------------------------
1351 subroutine bc_mxll_initialize_pml_simple(bc, space, gr, c_factor, dt)
1352 type(bc_mxll_t), intent(inout) :: bc
1353 class(space_t), intent(in) :: space
1354 type(grid_t), intent(in) :: gr
1355 real(real64), intent(in) :: c_factor
1356 real(real64), intent(in) :: dt
1357
1358 integer :: ip_in, ip, idir
1359 real(real64) :: ddv(1:space%dim), kappa, sigma, alpha
1361
1362 ! PML variables for all boundary points
1363 do ip_in = 1, bc%pml%points_number
1364 ip = bc%pml%points_map(ip_in)
1365 ddv(1:3) = abs(gr%x(ip, 1:3)) - bc%bc_bounds(1, 1:3)
1366 do idir = 1, space%dim
1367 if (ddv(idir) >= m_zero) then
1368 sigma = (ddv(idir)/bc%pml%width)**bc%pml%power * &
1369 (-(bc%pml%power + m_one)*p_c*c_factor*p_ep*log(bc%pml%refl_error)/(m_two * bc%pml%width))
1370 ! kappa and alpha could be set to different values, but these seem to be fine
1371 kappa = m_one
1372 ! non-zero values for alpha can be used to modify the low-frequency behavior
1373 alpha = m_zero
1374 bc%pml%b(ip_in, idir) = exp(-(alpha + sigma/kappa)/p_ep*dt)
1375 if (abs(sigma) < m_epsilon) then
1376 bc%pml%c(ip_in, idir) = m_zero
1377 else
1378 bc%pml%c(ip_in, idir) = m_one/(kappa + kappa**2*alpha/sigma) * &
1379 (bc%pml%b(ip_in, idir) - 1)
1380 end if
1381 else
1382 bc%pml%b(ip_in, idir) = m_zero
1383 bc%pml%c(ip_in, idir) = m_zero
1384 end if
1385 end do
1386 end do
1387 if (accel_is_enabled()) then
1388 call accel_write_buffer(bc%pml%buff_map, bc%pml%points_number, bc%pml%points_map)
1389 call accel_write_buffer(bc%pml%buff_b, bc%pml%points_number, space%dim, bc%pml%b)
1390 call accel_write_buffer(bc%pml%buff_c, bc%pml%points_number, space%dim, bc%pml%c)
1391 end if
1393 end subroutine bc_mxll_initialize_pml_simple
1394
1395 ! ---------------------------------------------------------
1396 subroutine bc_mxll_generate_mask(bc, mesh, bounds)
1397 type(bc_mxll_t), intent(inout) :: bc
1398 class(mesh_t), intent(in) :: mesh
1399 real(real64), intent(in) :: bounds(:,:)
1400
1401 integer :: ip, ip_in, idim, ip_in_max
1402 real(real64) :: ddv(3), tmp(3), width(3)
1403 real(real64), allocatable :: mask(:)
1404
1405 push_sub(bc_mxll_generate_mask)
1406
1407 call profiling_in('BC_MXLL_GENERATE_MASK')
1408
1409 ip_in_max = maxval(bc%mask_points_number(:))
1410
1411 safe_allocate(mask(1:mesh%np))
1412
1413 mask(:) = m_one
1414
1415 width(1:3) = bounds(2, 1:3) - bounds(1, 1:3)
1416 tmp(:) = m_zero
1417
1418 do ip = 1, mesh%np
1419 tmp = m_one
1420 mask(ip) = m_one
1421 ddv(1:3) = abs(mesh%x(ip, 1:3)) - bounds(1, 1:3)
1422 do idim = 1, mesh%box%dim
1423 if (ddv(idim) >= m_zero) then
1424 if (ddv(idim) <= width(idim)) then
1425 tmp(idim) = m_one - sin(ddv(idim) * m_pi / (m_two * (width(idim))))**2
1426 else
1427 tmp(idim) = m_one
1428 end if
1429 end if
1430 mask(ip) = mask(ip) * tmp(idim)
1431 end do
1432 end do
1433
1434 do idim = 1, mesh%box%dim
1435 do ip_in = 1, bc%mask_points_number(idim)
1436 ip = bc%mask_points_map(ip_in, idim)
1437 bc%mask(ip_in,idim) = mask(ip)
1438 end do
1439 end do
1440
1441 safe_deallocate_a(mask)
1442
1443 call profiling_out('BC_MXLL_GENERATE_MASK')
1444
1445 pop_sub(bc_mxll_generate_mask)
1447
1448 ! ---------------------------------------------------------
1449 subroutine bc_mxll_generate_medium(bc, space, gr, bounds, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
1450 type(bc_mxll_t), intent(inout) :: bc
1451 class(space_t), intent(in) :: space
1452 type(grid_t), intent(in) :: gr
1453 real(real64), intent(in) :: bounds(:,:)
1454 real(real64), intent(in) :: ep_factor
1455 real(real64), intent(in) :: mu_factor
1456 real(real64), intent(in) :: sigma_e_factor
1457 real(real64), intent(in) :: sigma_m_factor
1458
1459 integer :: ip, ipp, ip_in, ip_in_max, ip_bd, idim, point_info
1460 real(real64) :: dd, dd_min, dd_max, xx(3), xxp(3)
1461 real(real64), allocatable :: tmp(:), tmp_grad(:,:)
1462
1463 push_sub(bc_mxll_generate_medium)
1464
1465 call profiling_in('BC_MXLL_GENERATE_MEDIUM')
1466
1467 ip_in_max = max(bc%medium(1)%points_number, bc%medium(2)%points_number, bc%medium(3)%points_number)
1468 dd_max = max(2*gr%spacing(1), 2*gr%spacing(2), 2*gr%spacing(3))
1469
1470 do idim = 1, 3
1471 call single_medium_box_allocate(bc%medium(idim), ip_in_max)
1472 safe_allocate(tmp(gr%np_part))
1473 safe_allocate(tmp_grad(1:gr%np_part,1:space%dim))
1474 bc%medium(idim)%aux_ep = m_zero
1475 bc%medium(idim)%aux_mu = m_zero
1476 bc%medium(idim)%c = p_c
1477
1478 tmp = p_ep
1479 do ip = 1, gr%np_part
1480 call maxwell_box_point_info(bc, gr, ip, bounds, point_info)
1481 if ((point_info /= 0) .and. (abs(gr%x(ip, idim)) <= bounds(1, idim))) then
1482 xx(:) = gr%x(ip, :)
1483 dd_min = m_huge
1484 do ip_bd = 1, bc%medium(idim)%bdry_number
1485 ipp = bc%medium(idim)%bdry_map(ip_bd)
1486 xxp(:) = gr%x(ipp, :)
1487 dd = norm2(xx(1:3) - xxp(1:3))
1488 if (dd < dd_min) dd_min = dd
1489 end do
1490 tmp(ip) = p_ep * (m_one + ep_factor * m_one/(m_one + exp(-m_five/dd_max * (dd_min-m_two*dd_max))))
1491 end if
1492 end do
1493 call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
1494 do ip_in = 1, bc%medium(idim)%points_number
1495 ip = bc%medium(idim)%points_map(ip_in)
1496 bc%medium(idim)%aux_ep(ip_in, :) = &
1497 tmp_grad(ip, :)/(m_four*p_ep*ep_factor * m_one/(m_one + exp(-m_five/dd_max-dd)))
1498 end do
1499 end do
1500
1501 do idim = 1, 3
1502 tmp = p_mu
1503 do ip = 1, gr%np_part
1504 call maxwell_box_point_info(bc, gr, ip, bounds, point_info)
1505 if ((point_info == 1) .and. (abs(gr%x(ip, idim)) <= bounds(1, idim))) then
1506 xx(:) = gr%x(ip, :)
1507 dd_min = m_huge
1508 do ip_bd = 1, bc%medium(idim)%bdry_number
1509 ipp = bc%medium(idim)%bdry_map(ip_bd)
1510 xxp(:) = gr%x(ipp,:)
1511 dd = norm2(xx(1:3) - xxp(1:3))
1512 if (dd < dd_min) dd_min = dd
1513 end do
1514 tmp(ip) = p_mu * (m_one + mu_factor * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max))))
1515 end if
1516 end do
1517 call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
1518 do ip_in = 1, bc%medium(idim)%points_number
1519 ip = bc%medium(idim)%points_map(ip_in)
1520 bc%medium(idim)%aux_mu(ip_in, :) = &
1521 tmp_grad(ip, :)/(m_four * p_mu * mu_factor * m_one/(m_one + exp(-m_five/dd_max-dd)))
1522 end do
1523 end do
1524
1525 do idim = 1, 3
1526 do ip_in = 1, bc%medium(idim)%points_number
1527 ip = bc%medium(idim)%points_map(ip_in)
1528 xx(:) = gr%x(ip, :)
1529 dd_min = m_huge
1530 do ip_bd = 1, bc%medium(idim)%bdry_number
1531 ipp = bc%medium(idim)%bdry_map(ip_bd)
1532 xxp(:) = gr%x(ipp, :)
1533 dd = norm2(xx(1:3) - xxp(1:3))
1534 if (dd < dd_min) dd_min = dd
1535 end do
1536 bc%medium(idim)%ep(ip_in) = p_ep * (m_one + ep_factor &
1537 * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
1538 bc%medium(idim)%mu(ip_in) = p_mu * (m_one + mu_factor &
1539 * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
1540 bc%medium(idim)%sigma_e(ip_in) = (m_one + sigma_e_factor &
1541 * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
1542 bc%medium(idim)%sigma_m(ip_in) = (m_one + sigma_m_factor &
1543 * m_one/(m_one + exp(-m_five / dd_max * (dd_min - m_two * dd_max))))
1544 bc%medium(idim)%c(ip_in) = m_one / sqrt(bc%medium(idim)%ep(ip_in) * bc%medium(idim)%mu(ip_in))
1545 end do
1546 end do
1547
1548 safe_deallocate_a(tmp)
1549 safe_deallocate_a(tmp_grad)
1550
1551 call profiling_out('BC_MXLL_GENERATE_MEDIUM')
1552
1554 end subroutine bc_mxll_generate_medium
1555
1556 ! ---------------------------------------------------------
1557 subroutine maxwell_surfaces_init(mesh, st, bounds)
1558 class(mesh_t), intent(in) :: mesh
1559 type(states_mxll_t), intent(inout) :: st
1560 real(real64), intent(in) :: bounds(:,:)
1561
1562
1563 push_sub(maxwell_surfaces_init)
1564
1565 call profiling_in('MAXWELL_SURFACES_INIT')
1566
1567 ! y-z surface at -x boundary
1568 st%surface(1, 1)%spacing = m_half*(mesh%spacing(2) + mesh%spacing(3))
1569 st%surface(1, 1)%origin(:) = m_zero
1570 st%surface(1, 1)%origin(1) = -bounds(1, 1)
1571 st%surface(1, 1)%n(:) = m_zero
1572 st%surface(1, 1)%n(1) = -m_one
1573 st%surface(1, 1)%u(:) = m_zero
1574 st%surface(1, 1)%u(2) = -m_one
1575 st%surface(1, 1)%v(:) = m_zero
1576 st%surface(1, 1)%v(3) = m_one
1577 st%surface(1, 1)%nu = -int(bounds(1, 2)/mesh%spacing(2))
1578 st%surface(1, 1)%mu = int(bounds(1, 2)/mesh%spacing(2))
1579 st%surface(1, 1)%nv = -int(bounds(1, 3)/mesh%spacing(3))
1580 st%surface(1, 1)%mv = int(bounds(1, 3)/mesh%spacing(3))
1581
1582 ! y-z surface at +x boundary
1583 st%surface(2, 1)%spacing = m_half*(mesh%spacing(2) + mesh%spacing(3))
1584 st%surface(2, 1)%origin(:) = m_zero
1585 st%surface(2, 1)%origin(1) = bounds(1, 1)
1586 st%surface(2, 1)%n(:) = m_zero
1587 st%surface(2, 1)%n(1) = m_one
1588 st%surface(2, 1)%u(:) = m_zero
1589 st%surface(2, 1)%u(2) = m_one
1590 st%surface(2, 1)%v(:) = m_zero
1591 st%surface(2, 1)%v(3) = m_one
1592 st%surface(2, 1)%nu = -int(bounds(1, 2)/mesh%spacing(2))
1593 st%surface(2, 1)%mu = int(bounds(1, 2)/mesh%spacing(2))
1594 st%surface(2, 1)%nv = -int(bounds(1, 3)/mesh%spacing(3))
1595 st%surface(2, 1)%mv = int(bounds(1, 3)/mesh%spacing(3))
1596
1597 ! x-z surface at -y boundary
1598 st%surface(1, 2)%spacing = m_half*(mesh%spacing(1) + mesh%spacing(3))
1599 st%surface(1, 2)%origin(:) = m_zero
1600 st%surface(1, 2)%origin(2) = -bounds(1, 2)
1601 st%surface(1, 2)%n(:) = m_zero
1602 st%surface(1, 2)%n(2) = -m_one
1603 st%surface(1, 2)%u(:) = m_zero
1604 st%surface(1, 2)%u(1) = m_one
1605 st%surface(1, 2)%v(:) = m_zero
1606 st%surface(1, 2)%v(3) = m_one
1607 st%surface(1, 2)%nu = -int(bounds(1, 1)/mesh%spacing(1))
1608 st%surface(1, 2)%mu = int(bounds(1, 1)/mesh%spacing(1))
1609 st%surface(1, 2)%nv = -int(bounds(1, 3)/mesh%spacing(3))
1610 st%surface(1, 2)%mv = int(bounds(1, 3)/mesh%spacing(3))
1611
1612 ! x-z surface at +y boundary
1613 st%surface(2, 2)%spacing = m_half*(mesh%spacing(1) + mesh%spacing(3))
1614 st%surface(2, 2)%origin(:) = m_zero
1615 st%surface(2, 2)%origin(2) = bounds(1, 2)
1616 st%surface(2, 2)%n(:) = m_zero
1617 st%surface(2, 2)%n(2) = m_one
1618 st%surface(2, 2)%u(:) = m_zero
1619 st%surface(2, 2)%u(1) = m_one
1620 st%surface(2, 2)%v(:) = m_zero
1621 st%surface(2, 2)%v(3) = -m_one
1622 st%surface(2, 2)%nu = -int(bounds(1, 1)/mesh%spacing(1))
1623 st%surface(2, 2)%mu = int(bounds(1, 1)/mesh%spacing(1))
1624 st%surface(2, 2)%nv = -int(bounds(1, 3)/mesh%spacing(3))
1625 st%surface(2, 2)%mv = int(bounds(1, 3)/mesh%spacing(3))
1626
1627 ! x-y surface at -z boundary
1628 st%surface(1, 3)%spacing = m_half*(mesh%spacing(1) + mesh%spacing(2))
1629 st%surface(1, 3)%origin(:) = m_zero
1630 st%surface(1, 3)%origin(3) = -bounds(1, 3)
1631 st%surface(1, 3)%n(:) = m_zero
1632 st%surface(1, 3)%n(3) = -m_one
1633 st%surface(1, 3)%u(:) = m_zero
1634 st%surface(1, 3)%u(1) = m_one
1635 st%surface(1, 3)%v(:) = m_zero
1636 st%surface(1, 3)%v(2) = -m_one
1637 st%surface(1, 3)%nu = -int(bounds(1, 1)/mesh%spacing(1))
1638 st%surface(1, 3)%mu = int(bounds(1, 1)/mesh%spacing(1))
1639 st%surface(1, 3)%nv = -int(bounds(1, 2)/mesh%spacing(2))
1640 st%surface(1, 3)%mv = int(bounds(1, 2)/mesh%spacing(2))
1641
1642 ! x-y surface at +z boundary
1643 st%surface(2, 3)%spacing = m_half*(mesh%spacing(1) + mesh%spacing(2))
1644 st%surface(2, 3)%origin(:) = m_zero
1645 st%surface(2, 3)%origin(3) = bounds(1, 3)
1646 st%surface(2, 3)%n(:) = m_zero
1647 st%surface(2, 3)%n(3) = m_one
1648 st%surface(2, 3)%u(:) = m_zero
1649 st%surface(2, 3)%u(1) = m_one
1650 st%surface(2, 3)%v(:) = m_zero
1651 st%surface(2, 3)%v(2) = m_one
1652 st%surface(2, 3)%nu = -int(bounds(1, 1)/mesh%spacing(1))
1653 st%surface(2, 3)%mu = int(bounds(1, 1)/mesh%spacing(1))
1654 st%surface(2, 3)%nv = -int(bounds(1, 2)/mesh%spacing(2))
1655 st%surface(2, 3)%mv = int(bounds(1, 2)/mesh%spacing(2))
1656
1657 call profiling_out('MAXWELL_SURFACES_INIT')
1658
1659 pop_sub(maxwell_surfaces_init)
1660 end subroutine maxwell_surfaces_init
1661
1662 ! ---------------------------------------------------------
1663 subroutine maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
1664 type(bc_mxll_t), intent(inout) :: bc
1665 class(mesh_t), intent(in) :: mesh
1666 integer, intent(in) :: ip
1667 real(real64), intent(in) :: bounds(:,:)
1668 integer, intent(out) :: point_info
1669
1670 real(real64) :: rr, dd, xx(3), width(3)
1671
1672 point_info = 0
1673
1674 width(1:3) = bounds(2, 1:3) - bounds(1, 1:3)
1675 xx = m_zero
1676 xx(1:mesh%box%dim) = mesh%x(ip, 1:mesh%box%dim)
1677
1678 if (bc%ab_user_def) then
1679
1680 dd = bc%ab_ufn(ip) - bounds(1, 1)
1681 if (dd > m_zero) then
1682 if (bc%ab_ufn(ip) < bounds(2, 1)) then
1683 point_info = 1
1684 end if
1685 end if
1686
1687 else ! bc%ab_user_def == .false.
1688
1689 select type (box => mesh%box)
1690 type is (box_sphere_t)
1691 rr = norm2(xx - box%center)
1692 dd = rr - bounds(1, 1)
1693 if (dd > m_zero) then
1694 if (dd < width(1)) then
1695 point_info = 1
1696 end if
1697 end if
1698
1699 type is (box_parallelepiped_t)
1700 ! Limits of boundary region
1701 if (all(abs(xx(1:3) - box%center) <= bounds(2, 1:3))) then
1702 if (any(abs(xx(1:3) - box%center) > bounds(1, 1:3))) then
1703 point_info = 1
1704 else
1705 point_info = 0
1706 end if
1707 else
1708 point_info = -1
1709 end if
1710
1711 class default
1712 ! Other box shapes are not supported.
1713 assert(.false.)
1714 end select
1715 end if
1716
1717 end subroutine maxwell_box_point_info
1718
1719 ! ---------------------------------------------------------
1720 subroutine maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
1721 class(mesh_t), intent(in) :: mesh
1722 integer, intent(in) :: ip
1723 real(real64), intent(in) :: bounds(:,:)
1724 integer, intent(out) :: boundary_info
1725
1726 real(real64) :: xx(3)
1727
1728 boundary_info = 0
1729
1730 xx = m_zero
1731 xx(1:mesh%box%dim) = mesh%x(ip, 1:mesh%box%dim)
1732 if (is_close(abs(xx(1)), bounds(1, 1)) .and. (all(abs(xx(2:3)) <= bounds(1, 2:3))) .or. &
1733 is_close(abs(xx(2)), bounds(1, 2)) .and. (all(abs(xx(1:3:2)) <= bounds(1, 1:3:2))) .or. &
1734 is_close(abs(xx(3)), bounds(1, 3)) .and. (all(abs(xx(1:2)) <= bounds(1, 1:2)))) then
1735 boundary_info = 1
1736 else
1737 boundary_info = 0
1738 end if
1739
1740 end subroutine maxwell_boundary_point_info
1741
1742 ! ---------------------------------------------------------
1743 subroutine inner_and_outer_points_mapping(mesh, st, bounds)
1744 class(mesh_t), intent(in) :: mesh
1745 type(states_mxll_t), intent(inout) :: st
1746 real(real64), intent(in) :: bounds(:,:)
1747
1748 integer :: ip, ip_in, ip_bd, point_info
1749 real(real64) :: xx(mesh%box%dim)
1750
1752
1753 call profiling_in('IN_OUT_POINTS_MAP')
1754
1755 ! allocate inner and boundary points points map
1756 ip_in = 0
1757 ip_bd = 0
1758 do ip = 1, mesh%np
1759 xx = mesh%x(ip, :)
1760 if (all(abs(xx) <= bounds(2,1:mesh%box%dim))) then
1761 if (any(abs(xx) > bounds(1,1:mesh%box%dim))) then
1762 point_info = 1
1763 else
1764 point_info = 0
1765 end if
1766 else
1767 point_info = -1
1768 end if
1769 if (point_info == 0) then
1770 ip_in = ip_in + 1
1771 else
1772 ip_bd = ip_bd + 1
1773 end if
1774 end do
1775 st%inner_points_number = ip_in
1776 safe_allocate(st%inner_points_map(1:ip_in))
1777 safe_allocate(st%inner_points_mask(1:mesh%np))
1778 st%boundary_points_number = ip_bd
1779 safe_allocate(st%boundary_points_map(1:ip_bd))
1780 safe_allocate(st%boundary_points_mask(1:mesh%np))
1781 st%inner_points_mask = .false.
1782 st%boundary_points_mask = .false.
1783
1784 ! inner and boundary points mapping
1785 ip_in = 0
1786 ip_bd = 0
1787 do ip = 1, mesh%np
1788 xx = mesh%x(ip, :)
1789 if (all(abs(xx) <= bounds(2,1:mesh%box%dim))) then
1790 if (any(abs(xx) > bounds(1,1:mesh%box%dim))) then
1791 point_info = 1
1792 else
1793 point_info = 0
1794 end if
1795 else
1796 point_info = -1
1797 end if
1798 if (point_info == 0) then
1799 ip_in = ip_in + 1
1800 st%inner_points_map(ip_in) = ip
1801 st%inner_points_mask(ip) = .true.
1802 else
1803 ip_bd = ip_bd + 1
1804 st%boundary_points_map(ip_bd) = ip
1805 st%boundary_points_mask(ip) = .true.
1806 end if
1807 end do
1808
1809 if (accel_is_enabled()) then
1810 call accel_create_buffer(st%buff_inner_points_map, accel_mem_read_only, type_integer, st%inner_points_number)
1811 call accel_write_buffer(st%buff_inner_points_map, st%inner_points_number, st%inner_points_map)
1812 call accel_create_buffer(st%buff_boundary_points_map, accel_mem_read_only, type_integer, st%boundary_points_number)
1813 call accel_write_buffer(st%buff_boundary_points_map, st%boundary_points_number, st%boundary_points_map)
1814 end if
1816 call profiling_out('IN_OUT_POINTS_MAP')
1817
1819 end subroutine inner_and_outer_points_mapping
1820
1821 ! ---------------------------------------------------------
1822 subroutine surface_grid_points_mapping(mesh, st, bounds)
1823 class(mesh_t), intent(in) :: mesh
1824 type(states_mxll_t), intent(inout) :: st
1825 real(real64), intent(in) :: bounds(:,:)
1826
1827 integer :: ix, ix_max, iix, iy, iy_max, iiy, iz, iz_max, iiz, idx1, idx2, nn_max
1828 integer(int64) :: ip_global
1829 integer, allocatable :: nn(:,:,:,:)
1830 real(real64) :: rr(3), delta(3), vec(2), min_1(3), max_1(3), min_2(3), max_2(3)
1831
1833
1834 call profiling_in('SURF_GRID_POINTS_MAPPING')
1835
1836 st%surface_grid_rows_number(1) = 3
1837 ix_max = st%surface_grid_rows_number(1)
1838 st%surface_grid_rows_number(2) = 3
1839 iy_max = st%surface_grid_rows_number(2)
1840 st%surface_grid_rows_number(3) = 3
1841 iz_max = st%surface_grid_rows_number(3)
1842
1843 delta(1) = m_two * abs(bounds(1,1)) / real(ix_max, real64)
1844 delta(2) = m_two * abs(bounds(1,2)) / real(iy_max, real64)
1845 delta(3) = m_two * abs(bounds(1,3)) / real(iz_max, real64)
1846
1847 st%surface_grid_element(1) = delta(2) * delta(3)
1848 st%surface_grid_element(2) = delta(1) * delta(3)
1849 st%surface_grid_element(3) = delta(1) * delta(2)
1850
1851 safe_allocate(nn(1:2, 1:3, 1:3, 1:3))
1852
1853 st%surface_grid_center(1, 1, :, :) = -int(bounds(1,1))
1854 do iy = 1, iy_max
1855 do iz = 1, iz_max
1856 rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
1857 rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
1858 st%surface_grid_center(1, 2, iy, iz) = int(rr(2))
1859 st%surface_grid_center(1, 3, iy, iz) = int(rr(3))
1860 end do
1861 end do
1862 st%surface_grid_center(2, 1, :, :) = int(bounds(1,1))
1863 do iy = 1, iy_max
1864 do iz = 1, iz_max
1865 rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
1866 rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
1867 st%surface_grid_center(2, 2, iy, iz) = int(rr(2))
1868 st%surface_grid_center(2, 3, iy, iz) = int(rr(3))
1869 end do
1870 end do
1871
1872 st%surface_grid_center(1, 2, :, :) = -int(bounds(1,2))
1873 do ix = 1, ix_max
1874 do iz = 1, iz_max
1875 rr(1) = -bounds(1,1) + delta(1)/m_two + (ix-1) * delta(1)
1876 rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
1877 st%surface_grid_center(1, 1, ix, iz) = int(rr(1))
1878 st%surface_grid_center(1, 3, ix, iz) = int(rr(3))
1879 end do
1880 end do
1881 st%surface_grid_center(2, 2, :, :) = int(bounds(1,2))
1882 do ix = 1, ix_max
1883 do iz = 1, iz_max
1884 rr(1) = -bounds(1,2) + delta(1)/m_two + (ix-1) * delta(1)
1885 rr(3) = -bounds(1,3) + delta(3)/m_two + (iz-1) * delta(3)
1886 st%surface_grid_center(2, 1, ix, iz) = int(rr(1))
1887 st%surface_grid_center(2, 3, ix, iz) = int(rr(3))
1888 end do
1889 end do
1890
1891 st%surface_grid_center(1, 3, :, :) = -int(bounds(1,3))
1892 do ix = 1, ix_max
1893 do iy = 1, iy_max
1894 rr(1) = -bounds(1,1) + delta(1)/m_two + (ix-1) * delta(1)
1895 rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
1896 st%surface_grid_center(1, 1, ix, iy) = int(rr(1))
1897 st%surface_grid_center(1, 2, ix, iy) = int(rr(2))
1898 end do
1899 end do
1900 st%surface_grid_center(2, 3, :, :) = int(bounds(1,3))
1901 do ix = 1, ix_max
1902 do iy = 1, iy_max
1903 rr(1) = -bounds(1,2) + delta(1)/m_two + (ix-1) * delta(1)
1904 rr(2) = -bounds(1,2) + delta(2)/m_two + (iy-1) * delta(2)
1905 st%surface_grid_center(2, 1, ix, iy) = int(rr(1))
1906 st%surface_grid_center(2, 2, ix, iy) = int(rr(2))
1907 end do
1908 end do
1909
1910 st%surface_grid_points_number(:,:,:) = 0
1911
1912 nn_max = 0
1913
1914 do iy = 1, iy_max
1915 do iz = 1, iz_max
1916 min_1(iy) = -bounds(1,2) + (iy-1) * delta(2)
1917 max_1(iy) = -bounds(1,2) + iy * delta(2)
1918 min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
1919 max_2(iz) = -bounds(1,3) + iz * delta(3)
1920 end do
1921 end do
1922 do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
1923 do iiz = mesh%idx%nr(1,3), mesh%idx%nr(2,3)
1924 vec(1) = iiy * mesh%spacing(2)
1925 vec(2) = iiz * mesh%spacing(3)
1926 call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
1927 if (idx1 /= 0 .and. idx2 /= 0) then
1928 st%surface_grid_points_number(1, idx1, idx2) = st%surface_grid_points_number(1, idx1, idx2) + 1
1929 if (nn_max < st%surface_grid_points_number(1, idx1, idx2)) then
1930 nn_max = st%surface_grid_points_number(1, idx1, idx2)
1931 end if
1932 end if
1933 end do
1934 end do
1935
1936 do ix = 1, ix_max
1937 do iz = 1, iz_max
1938 min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
1939 max_1(ix) = -bounds(1,1) + ix * delta(1)
1940 min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
1941 max_2(iz) = -bounds(1,3) + iz * delta(3)
1942 end do
1943 end do
1944 do iix = mesh%idx%nr(1, 1), mesh%idx%nr(2, 1)
1945 do iiz = mesh%idx%nr(1, 3), mesh%idx%nr(2, 3)
1946 vec(1) = iix * mesh%spacing(1)
1947 vec(2) = iiz * mesh%spacing(3)
1948 call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
1949 if (idx1 /= 0 .and. idx2 /= 0) then
1950 st%surface_grid_points_number(2,idx1,idx2) = st%surface_grid_points_number(2,idx1,idx2)+1
1951 if (nn_max < st%surface_grid_points_number(2, idx1, idx2)) then
1952 nn_max = st%surface_grid_points_number(2, idx1, idx2)
1953 end if
1954 end if
1955 end do
1956 end do
1957
1958 do ix = 1, ix_max
1959 do iy = 1, iy_max
1960 min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
1961 max_1(ix) = -bounds(1,1) + ix * delta(1)
1962 min_2(iy) = -bounds(1,2) + (iy-1) * delta(2)
1963 max_2(iy) = -bounds(1,2) + iy * delta(2)
1964 end do
1965 end do
1966 do iix = mesh%idx%nr(1,1), mesh%idx%nr(2,1)
1967 do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
1968 vec(1) = iix * mesh%spacing(1)
1969 vec(2) = iiy * mesh%spacing(2)
1970 call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
1971 if (idx1 /= 0 .and. idx2 /= 0) then
1972 st%surface_grid_points_number(3, idx1, idx2) = st%surface_grid_points_number(3, idx1, idx2) + 1
1973 if (nn_max < st%surface_grid_points_number(3, idx1, idx2)) then
1974 nn_max = st%surface_grid_points_number(3, idx1, idx2)
1975 end if
1976 end if
1977 end do
1978 end do
1979
1980 ! originally there were three allocated of the same pointer here
1981 safe_allocate(st%surface_grid_points_map(1:2, 1:st%dim, 1:ix_max, 1:iy_max, 1:nn_max))
1982
1983 nn(:,:,:,:) = 0
1984
1985 do iy = 1, iy_max
1986 do iz = 1, iz_max
1987 min_1(iy) = -bounds(1,2) + (iy-1) * delta(2)
1988 max_1(iy) = -bounds(1,2) + iy * delta(2)
1989 min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
1990 max_2(iz) = -bounds(1,3) + iz * delta(3)
1991 end do
1992 end do
1993 do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
1994 do iiz = mesh%idx%nr(1,3), mesh%idx%nr(2,3)
1995 vec(1) = iiy * mesh%spacing(2)
1996 vec(2) = iiz * mesh%spacing(3)
1997 call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
1998 if (idx1 /= 0 .and. idx2 /= 0) then
1999 nn(1, 1, idx1, idx2) = nn(1, 1, idx1, idx2) + 1
2000 rr(1) = -bounds(1, 1)
2001 rr(2) = iiy * mesh%spacing(2)
2002 rr(3) = iiz * mesh%spacing(3)
2003 iix = int(-bounds(1,1)/mesh%spacing(1))
2004 ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
2005 st%surface_grid_points_map(1, 1, idx1, idx2, nn(1, 1, idx1, idx2)) = ip_global
2006 nn(2, 1, idx1, idx2) = nn(2, 1, idx1, idx2) + 1
2007 rr(1) = bounds(1,1)
2008 rr(2) = iiy * mesh%spacing(2)
2009 rr(3) = iiz * mesh%spacing(3)
2010 iix = int(bounds(1,1)/mesh%spacing(1))
2011 ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
2012 st%surface_grid_points_map(2, 1, idx1, idx2, nn(2, 1, idx1, idx2)) = ip_global
2013 end if
2014 end do
2015 end do
2016
2017 do ix = 1, ix_max
2018 do iz = 1, iz_max
2019 min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
2020 max_1(ix) = -bounds(1,1) + ix * delta(1)
2021 min_2(iz) = -bounds(1,3) + (iz-1) * delta(3)
2022 max_2(iz) = -bounds(1,3) + iz * delta(3)
2023 end do
2024 end do
2025 do iix = mesh%idx%nr(1,1), mesh%idx%nr(2,1)
2026 do iiz = mesh%idx%nr(1,3), mesh%idx%nr(2,3)
2027 vec(1) = iix * mesh%spacing(1)
2028 vec(2) = iiz * mesh%spacing(3)
2029 call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
2030 if (idx1 /= 0 .and. idx2 /= 0) then
2031 nn(1, 2, idx1, idx2) = nn(1, 2, idx1, idx2) + 1
2032 rr(1) = iix * mesh%spacing(1)
2033 rr(2) = -bounds(1, 2)
2034 rr(3) = iiz * mesh%spacing(3)
2035 iiy = int(-bounds(1,2)/mesh%spacing(2))
2036 ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
2037 st%surface_grid_points_map(1, 2, idx1, idx2, nn(1, 2, idx1, idx2)) = ip_global
2038 nn(2, 2, idx1, idx2) = nn(2, 2, idx1, idx2) + 1
2039 rr(1) = iix * mesh%spacing(1)
2040 rr(2) = bounds(1,2)
2041 rr(3) = iiz * mesh%spacing(3)
2042 iiy = int(bounds(1,2)/mesh%spacing(2))
2043 ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
2044 st%surface_grid_points_map(2, 2, idx1, idx2, nn(2, 2, idx1, idx2)) = ip_global
2045 end if
2046 end do
2047 end do
2048
2049 do ix = 1, ix_max
2050 do iy = 1, iy_max
2051 min_1(ix) = -bounds(1,1) + (ix-1) * delta(1)
2052 max_1(ix) = -bounds(1,1) + ix * delta(1)
2053 min_2(iy) = -bounds(1,2) + (iy-1) * delta(2)
2054 max_2(iy) = -bounds(1,2) + iy * delta(2)
2055 end do
2056 end do
2057 do iix = mesh%idx%nr(1,1), mesh%idx%nr(2,1)
2058 do iiy = mesh%idx%nr(1,2), mesh%idx%nr(2,2)
2059 vec(1) = iix * mesh%spacing(1)
2060 vec(2) = iiy * mesh%spacing(2)
2061 call get_surface_indices(vec, min_1, max_1, min_2, max_2, idx1, idx2)
2062 if (idx1 /= 0 .and. idx2 /= 0) then
2063 nn(1, 3, idx1, idx2) = nn(1, 3, idx1, idx2) + 1
2064 rr(1) = iix * mesh%spacing(1)
2065 rr(2) = iiy * mesh%spacing(2)
2066 rr(3) = -bounds(1,3)
2067 iiz = int(-bounds(1,3)/mesh%spacing(3))
2068 ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
2069 st%surface_grid_points_map(1, 3, idx1, idx2, nn(1, 3, idx1, idx2)) = ip_global
2070 nn(2, 3, idx1, idx2) = nn(2, 3, idx1, idx2) + 1
2071 rr(1) = iix * mesh%spacing(1)
2072 rr(2) = iiy * mesh%spacing(2)
2073 rr(3) = bounds(1,3)
2074 iiz = int(bounds(1,3)/mesh%spacing(3))
2075 ip_global = mesh_global_index_from_coords(mesh, [iix, iiy, iiz])
2076 st%surface_grid_points_map(2, 3, idx1, idx2, nn(2, 3, idx1, idx2)) = ip_global
2077 end if
2078 end do
2079 end do
2080
2081 safe_deallocate_a(nn)
2082
2083 call profiling_out('SURF_GRID_POINTS_MAPPING')
2084
2086 contains
2087
2088 subroutine get_surface_indices(vec, min_1, max_1, min_2, max_2, index_1, index_2)
2089 real(real64), intent(in) :: vec(:)
2090 real(real64), intent(in) :: min_1(:)
2091 real(real64), intent(in) :: max_1(:)
2092 real(real64), intent(in) :: min_2(:)
2093 real(real64), intent(in) :: max_2(:)
2094 integer, intent(out) :: index_1
2095 integer, intent(out) :: index_2
2096
2097 if (vec(1) >= min_1(1) .and. vec(1) <= max_1(1) .and. vec(2) >= min_2(1) .and. vec(2) <= max_2(1)) then
2098 index_1 = 1
2099 index_2 = 1
2100 else if (vec(1) >= min_1(2) .and. vec(1) <= max_1(2) .and. vec(2) >= min_2(1) .and. vec(2) <= max_2(1)) then
2101 index_1 = 2
2102 index_2 = 1
2103 else if (vec(1) >= min_1(3) .and. vec(1) <= max_1(3) .and. vec(2) >= min_2(1) .and. vec(2) <= max_2(1)) then
2104 index_1 = 3
2105 index_2 = 1
2106 else if (vec(1) >= min_1(1) .and. vec(1) <= max_1(1) .and. vec(2) >= min_2(2) .and. vec(2) <= max_2(2)) then
2107 index_1 = 1
2108 index_2 = 2
2109 else if (vec(1) >= min_1(2) .and. vec(1) <= max_1(2) .and. vec(2) >= min_2(2) .and. vec(2) <= max_2(2)) then
2110 index_1 = 2
2111 index_2 = 2
2112 else if (vec(1) >= min_1(3) .and. vec(1) <= max_1(3) .and. vec(2) >= min_2(2) .and. vec(2) <= max_2(2)) then
2113 index_1 = 3
2114 index_2 = 2
2115 else if (vec(1) >= min_1(1) .and. vec(1) <= max_1(1) .and. vec(2) >= min_2(3) .and. vec(2) <= max_2(3)) then
2116 index_1 = 1
2117 index_2 = 3
2118 else if (vec(1) >= min_1(2) .and. vec(1) <= max_1(2) .and. vec(2) >= min_2(3) .and. vec(2) <= max_2(3)) then
2119 index_1 = 2
2120 index_2 = 3
2121 else if (vec(1) >= min_1(3) .and. vec(1) <= max_1(3) .and. vec(2) >= min_2(3) .and. vec(2) <= max_2(3)) then
2122 index_1 = 3
2123 index_2 = 3
2124 else
2125 index_1 = 0
2126 index_2 = 0
2127 end if
2128
2129 end subroutine get_surface_indices
2130
2131 end subroutine surface_grid_points_mapping
2132
2134
2135!! Local Variables:
2136!! mode: f90
2137!! coding: utf-8
2138!! End:
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
subroutine get_surface_indices(vec, min_1, max_1, min_2, max_2, index_1, index_2)
subroutine get_medium_io_function(medium_func, bc, io_func, idim)
subroutine get_mask_io_function(mask_func, bc, io_func, idim)
subroutine write_files(filename, tmp)
subroutine get_pml_io_function(pml_func, bc, io_func)
integer, parameter, public accel_mem_read_write
Definition: accel.F90:196
subroutine, public accel_release_buffer(this, async)
Definition: accel.F90:919
pure logical function, public accel_is_enabled()
Definition: accel.F90:419
integer, parameter, public accel_mem_read_only
Definition: accel.F90:196
type(debug_t), save, public debug
Definition: debug.F90:158
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine, public external_waves_end(external_waves)
subroutine, public external_waves_init(external_waves, namespace)
Here, plane wave is evaluated from analytical formulae on grid.
real(real64), parameter, public m_two
Definition: global.F90:193
real(real64), parameter, public m_huge
Definition: global.F90:209
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
real(real64), parameter, public p_mu
Definition: global.F90:234
real(real64), parameter, public p_ep
Definition: global.F90:233
complex(real64), parameter, public m_z0
Definition: global.F90:201
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:229
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_five
Definition: global.F90:196
This module implements the underlying real-space grid.
Definition: grid.F90:119
This module implements the index, used for the mesh points.
Definition: index.F90:124
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
integer(int64) function, public io_function_fill_how(where)
Use this function to quickly plot functions for debugging purposes: call dio_function_output(io_funct...
Definition: io.F90:116
subroutine, public single_medium_box_allocate(medium_box, n_points)
Allocation of medium_box components.
subroutine, public single_medium_box_end(medium_box)
Deallocation of medium_box components.
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
integer, parameter, public mxll_ab_mask_zero
subroutine maxwell_zero_points_mapping(bc, mesh, bounds)
integer, parameter, public mxll_bc_mirror_pmc
integer, parameter mxll_plane_waves_positive_side
subroutine maxwell_boundary_point_info(mesh, ip, bounds, boundary_info)
subroutine, public bc_mxll_init(bc, namespace, space, gr, st)
subroutine bc_mxll_generate_medium(bc, space, gr, bounds, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
subroutine maxwell_surfaces_init(mesh, st, bounds)
subroutine maxwell_box_point_info(bc, mesh, ip, bounds, point_info)
integer, parameter, public mxll_bc_mirror_pec
integer, parameter, public mxll_ab_mask
subroutine bc_mxll_medium_init(gr, namespace, bounds, idim, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
subroutine, public surface_grid_points_mapping(mesh, st, bounds)
integer, parameter, public mxll_bc_plane_waves
subroutine, public bc_mxll_initialize_pml_simple(bc, space, gr, c_factor, dt)
integer, parameter, public mxll_bc_medium
subroutine maxwell_constant_points_mapping(bc, mesh, bounds)
subroutine bc_mxll_generate_pml(bc, space)
subroutine bc_mxll_write_info(bc, mesh, namespace, space)
subroutine, public bc_mxll_generate_pml_parameters(bc, space, gr, c_factor, dt)
subroutine bc_mxll_generate_mask(bc, mesh, bounds)
subroutine maxwell_plane_waves_points_mapping(bc, mesh, bounds)
subroutine, public inner_and_outer_points_mapping(mesh, st, bounds)
integer, parameter, public mxll_bc_constant
subroutine bc_mxll_pml_init(bc, gr, namespace, ab_bounds, idim)
subroutine maxwell_medium_points_mapping(bc, mesh, bounds)
subroutine bc_mxll_ab_bounds_init(bc, gr, namespace, ab_bounds, idim)
subroutine, public bc_mxll_end(bc)
integer, parameter, public mxll_ab_cpml
subroutine maxwell_pml_points_mapping(bc, mesh, bounds)
subroutine maxwell_mask_points_mapping(bc, mesh, bounds)
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
integer(int64) function, public mesh_global_index_from_coords(mesh, ix)
This function returns the true global index of the point for a given vector of integer coordinates.
Definition: mesh.F90:916
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:898
character(len=512), private msg
Definition: messages.F90:167
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
Some general things and nomenclature:
Definition: par_vec.F90:173
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:455
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
type(type_t), public type_float
Definition: types.F90:135
type(type_t), public type_cmplx
Definition: types.F90:136
type(type_t), public type_integer
Definition: types.F90:137
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:225
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing
type(unit_t), public unit_one
some special units required for particular quantities
Class implementing a parallelepiped box. Currently this is restricted to a rectangular cuboid (all th...
Class implementing a spherical box.
Definition: box_sphere.F90:134
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
int true(void)