67 integer :: points_number
68 integer,
allocatable :: points_map(:)
69 integer,
allocatable :: points_map_inv(:)
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.
87 type(accel_mem_t) :: buff_a, buff_b, buff_c, buff_conv_plus, buff_conv_minus, buff_map, buff_conv_plus_old
92 integer :: bc_ab_type(3)
93 real(real64) :: bc_bounds(2, 3)
94 logical :: ab_user_def
95 real(real64),
allocatable :: ab_ufn(:)
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(:,:)
104 type(single_medium_box_t) :: medium(3)
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
111 integer :: mirror_points_number(3)
112 integer,
allocatable :: mirror_points_map(:,:)
114 logical :: do_plane_waves = .false.
115 type(external_waves_t) :: plane_wave
116 logical :: plane_waves_dims(1:3) = .false.
118 real(real64) :: zero_width
119 integer :: zero_points_number(3)
120 integer,
allocatable :: zero_points_map(:,:)
121 real(real64),
allocatable :: zero(:,:)
124 integer,
public,
parameter :: &
132 integer,
parameter :: &
133 MXLL_PLANE_WAVES_NEGATIVE_SIDE = -1, &
136 integer,
public,
parameter :: &
137 MXLL_AB_NOT_ABSORBING = 0, &
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
152 integer :: idim, nlines, icol, ncols, ab_shape_dim
153 real(real64) :: bounds(2, space%dim), ab_bounds(2, space%dim)
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
166 plane_waves_set = .false.
167 bc%bc_type(:) = mxll_bc_zero
197 if (
parse_block(namespace,
'MaxwellBoundaryConditions', blk) == 0)
then
202 if (nlines /= 1)
then
216 select case (bc%bc_type(icol))
222 string =
'PEC Mirror'
224 string =
'PMC Mirror'
226 string =
'Plane waves'
228 string =
'Medium boundary'
230 write(
message(1),
'(a)')
'Unknown Maxwell boundary condition'
233 write(
message(1),
'(a,I1,a,a)')
'Maxwell boundary condition in direction ', icol,
': ', trim(string)
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!'
242 plane_waves_check = .false.
243 ab_mask_check = .false.
244 ab_pml_check = .false.
245 constant_check = .false.
248 bc%ab_user_def = .false.
249 bc%bc_ab_type(:) = mxll_ab_not_absorbing
273 if (
parse_block(namespace,
'MaxwellAbsorbingBoundaries', blk) == 0)
then
276 if (nlines /= 1)
then
277 message(1) =
'MaxwellAbsorbingBounaries has to consist of one line!'
283 message(1) =
'MaxwellAbsorbingBoundaries has to consist of three columns!'
298 if (bc%bc_type(idim) == mxll_bc_zero) zero_check = .
true.
301 if (ab_mask_check .or. ab_pml_check)
then
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."
312 bounds(1, :) = (gr%idx%nr(2, :) - gr%idx%enlarge(:))*gr%spacing(:)
313 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
315 select case (bc%bc_type(idim))
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)
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)
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.
335 call bc_mxll_medium_init(gr, namespace, bounds, idim, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
341 select type (box => gr%box)
344 if (space%is_periodic())
then
345 message(1) =
"Sphere box shape can only work for non-periodic systems"
349 if (.not. box%axes%orthogonal)
then
350 message(1) =
"Maxwell propagation does not work for non-orthogonal grids."
354 ab_shape_dim = space%dim
355 ab_bounds(1, idim) = bounds(1, idim)
356 ab_bounds(2, idim) = bounds(1, idim)
358 message(1) =
"Box shape for Maxwell propagation not supported yet"
362 if (bc%bc_ab_type(idim) /= mxll_ab_not_absorbing)
then
366 select case (bc%bc_ab_type(idim))
369 bc%zero_width = bc%ab_width
373 bc%mask_width = bc%ab_width
379 message(1) =
"Absorbing boundary type not implemented for Maxwell propagation"
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)
391 bc%bc_bounds(:, idim) = bounds(:, idim)
394 select type (box => gr%box)
396 select case (bc%bc_ab_type(idim))
409 string = trim(ab_type_str)//
" Absorbing Boundary"
410 write(string,
'(2a, f10.3,3a)') trim(string),
" in "//dims(idim)//
' direction spans from: ', &
413 write(
message(1),
'(a)') trim(string)
416 write(string,
'(a,f10.3,3a)') trim(string),&
420 write(
message(2),
'(a)') trim(string)
426 write(
message(1),
'(a,es10.3,3a)') &
429 write(
message(2),
'(a,es10.3,3a)') &
442 if (ab_mask_check)
then
447 if (ab_pml_check)
then
452 if (constant_check)
then
453 bounds(1, :) = gr%idx%nr(2, :) * gr%spacing(:)
454 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
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)
465 if (plane_waves_check)
then
466 bounds(1, :) = gr%idx%nr(2, :) * gr%spacing(:)
467 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
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)
480 bounds(1, :) = gr%idx%nr(2, :) * gr%spacing(:)
481 bounds(2, :) = gr%idx%nr(2, :) * gr%spacing(:)
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)
491 if (ab_mask_check)
then
495 if (ab_pml_check)
then
503 if (ab_mask_check .or. ab_pml_check)
then
520 safe_deallocate_a(bc%ab_ufn)
522 safe_deallocate_a(bc%mask_points_map)
523 safe_deallocate_a(bc%mask)
530 safe_deallocate_a(bc%constant_points_map)
531 safe_deallocate_a(bc%constant_rs_state)
536 safe_deallocate_a(bc%mirror_points_map)
540 safe_deallocate_a(bc%zero_points_map)
541 safe_deallocate_a(bc%zero)
548 type(
pml_t),
intent(inout) :: pml
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)
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
594 real(real64) :: width
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)
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
661 real(real64) :: default_width
673 default_width =
m_two * gr%der%order * maxval(gr%spacing(1:3))
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
683 ab_bounds(1, idim) = ab_bounds(2, idim) - bc%ab_width
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
699 bc%pml%width = bc%ab_width
718 call parse_variable(namespace,
'MaxwellABPMLReflectionError', 1.0e-16_real64, bc%pml%refl_error,
unit_one)
726 class(
mesh_t),
intent(in) :: mesh
727 type(namespace_t),
intent(in) :: namespace
728 class(
space_t),
intent(in) :: space
730 integer :: err, idim, idim2
731 real(real64),
allocatable :: tmp(:)
732 logical :: mask_check, pml_check, medium_check
733 character(1) :: dim_label(3)
741 medium_check = .false.
743 dim_label = (/
'x',
'y',
'z'/)
753 medium_check = .
true.
758 safe_allocate(tmp(1:mesh%np))
764 safe_deallocate_a(tmp)
765 else if (pml_check)
then
766 safe_allocate(tmp(1:mesh%np))
771 call write_files(
"maxwell_sigma_e-"//dim_label(idim), tmp)
776 call write_files(
"maxwell_sigma_m-"//dim_label(idim), tmp)
781 call write_files(
"maxwell_sigma_pml_a_e-"//dim_label(idim), tmp)
783 safe_deallocate_a(tmp)
786 if (medium_check)
then
787 safe_allocate(tmp(1:mesh%np))
792 call write_files(
"maxwell_ep"//dim_label(idim), tmp)
796 call write_files(
"maxwell_mu"//dim_label(idim), tmp)
800 call write_files(
"maxwell_c"//dim_label(idim), tmp)
805 call write_files(
"maxwell_aux_ep-"//dim_label(idim)//
"-"//dim_label(idim2), tmp)
810 call write_files(
"maxwell_aux_mu-"//dim_label(idim)//
"-"//dim_label(idim2), tmp)
814 safe_deallocate_a(tmp)
823 real(real64),
intent(in) :: pml_func(:)
825 real(real64),
intent(inout) :: io_func(:)
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)
837 real(real64),
intent(in) :: mask_func(:,:)
839 real(real64),
intent(inout) :: io_func(:)
840 integer,
intent(in) :: idim
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)
852 real(real64),
intent(in) :: medium_func(:)
854 real(real64),
intent(inout) :: io_func(:)
855 integer,
intent(in) :: idim
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)
867 character(len=*),
intent(in) :: filename
868 real(real64),
intent(in) :: tmp(:)
891 class(
mesh_t),
intent(in) :: mesh
892 real(real64),
intent(in) :: bounds(:,:)
894 integer :: ip, ip_in, ip_in_max, point_info, idim
906 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
910 if (ip_in > ip_in_max) ip_in_max = ip_in
911 bc%mask_points_number(idim) = ip_in
914 safe_allocate(bc%mask(1:ip_in_max, 1:idim))
915 safe_allocate(bc%mask_points_map(1:ip_in_max, 1:idim))
923 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
925 bc%mask_points_map(ip_in, idim) = ip
938 class(
mesh_t),
intent(in) :: mesh
939 real(real64),
intent(in) :: bounds(:,:)
941 integer :: ip, ip_in, point_info
951 if (point_info == 1)
then
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
962 if (point_info == 1)
then
964 bc%pml%points_map(ip_in) = ip
965 bc%pml%points_map_inv(ip) = ip_in
977 class(
mesh_t),
intent(in) :: mesh
978 real(real64),
intent(in) :: bounds(:,:)
980 integer :: ip, ip_in, point_info
990 if (point_info == 1)
then
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))
1002 if (point_info == 1)
then
1004 bc%constant_points_map(ip_in) = ip
1010 call accel_write_buffer(bc%buff_constant_points_map, bc%constant_points_number, bc%constant_points_map)
1021 class(
mesh_t),
intent(in) :: mesh
1022 real(real64),
intent(in) :: bounds(:,:)
1024 integer :: ip, ip_in, point_info
1036 if (point_info == 1)
then
1040 bc%plane_wave%points_number = ip_in
1041 safe_allocate(bc%plane_wave%points_map(1:ip_in))
1047 if (point_info == 1)
then
1049 bc%plane_wave%points_map(ip_in) = ip
1055 call accel_write_buffer(bc%plane_wave%buff_map, bc%plane_wave%points_number, bc%plane_wave%points_map)
1066 class(
mesh_t),
intent(in) :: mesh
1067 real(real64),
intent(in) :: bounds(:,:)
1069 integer :: ip, ip_in, ip_in_max, point_info, idim
1077 if (bc%bc_type(idim) == mxll_bc_zero)
then
1082 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
1086 if (ip_in > ip_in_max) ip_in_max = ip_in
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))
1094 if (bc%bc_type(idim) == mxll_bc_zero)
then
1099 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
1101 bc%zero_points_map(ip_in, idim) = ip
1115 class(
mesh_t),
intent(in) :: mesh
1116 real(real64),
intent(in) :: bounds(:,:)
1118 integer :: ip, ip_in, ip_in_max, ip_bd, ip_bd_max, point_info, boundary_info, idim
1134 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
1137 if ((boundary_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
1141 bc%medium(idim)%points_number = ip_in
1142 bc%medium(idim)%bdry_number = ip_bd
1146 safe_allocate(bc%medium(idim)%points_map(1:ip_in_max))
1147 safe_allocate(bc%medium(idim)%bdry_map(1:ip_bd_max))
1157 if ((point_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
1159 bc%medium(idim)%points_map(ip_in) = ip
1161 if ((boundary_info == 1) .and. (abs(mesh%x(ip, idim)) >= bounds(1, idim)))
then
1163 bc%medium(idim)%bdry_map(ip_bd) = ip
1177 class(
space_t),
intent(in) :: space
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))
1198 bc%pml%kappa =
m_one
1205 bc%pml%conv_plus =
m_z0
1206 bc%pml%conv_minus =
m_z0
1207 bc%pml%conv_plus_old =
m_z0
1216 int(bc%pml%points_number, int64)*space%dim*space%dim)
1218 int(bc%pml%points_number, int64)*space%dim*space%dim)
1220 int(bc%pml%points_number, int64)*space%dim*space%dim)
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)
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
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(:,:)
1251 safe_allocate(tmp(gr%np_part))
1252 safe_allocate(tmp_grad(gr%np, 1:space%dim))
1258 width(1:3) = (/ bc%pml%width, bc%pml%width, bc%pml%width /)
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
1269 ss_max = -(bc%pml%power +
m_one)*
p_c*c_factor*
p_ep*
log(bc%pml%refl_error)/(
m_two * bc%pml%width)
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)
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
1298 do idim = 1, space%dim
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)
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))
1312 do idim = 1, space%dim
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)
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))
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)
1332 safe_deallocate_a(tmp)
1333 safe_deallocate_a(tmp_grad)
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)
1344 bc%pml%parameters_initialized = .
true.
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
1358 integer :: ip_in, ip, idir
1359 real(real64) :: ddv(1:space%dim), kappa, sigma, alpha
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 * &
1374 bc%pml%b(ip_in, idir) =
exp(-(alpha + sigma/kappa)/
p_ep*dt)
1376 bc%pml%c(ip_in, idir) =
m_zero
1378 bc%pml%c(ip_in, idir) =
m_one/(kappa + kappa**2*alpha/sigma) * &
1379 (bc%pml%b(ip_in, idir) - 1)
1382 bc%pml%b(ip_in, idir) =
m_zero
1383 bc%pml%c(ip_in, idir) =
m_zero
1398 class(
mesh_t),
intent(in) :: mesh
1399 real(real64),
intent(in) :: bounds(:,:)
1401 integer :: ip, ip_in, idim, ip_in_max
1402 real(real64) :: ddv(3), tmp(3), width(3)
1403 real(real64),
allocatable :: mask(:)
1409 ip_in_max = maxval(bc%mask_points_number(:))
1411 safe_allocate(mask(1:mesh%np))
1415 width(1:3) = bounds(2, 1:3) - bounds(1, 1:3)
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
1430 mask(ip) = mask(ip) * tmp(idim)
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)
1441 safe_deallocate_a(mask)
1449 subroutine bc_mxll_generate_medium(bc, space, gr, bounds, ep_factor, mu_factor, sigma_e_factor, sigma_m_factor)
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
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(:,:)
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))
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
1479 do ip = 1, gr%np_part
1481 if ((point_info /= 0) .and. (abs(gr%x(ip, idim)) <= bounds(1, idim)))
then
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
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, :) = &
1503 do ip = 1, gr%np_part
1505 if ((point_info == 1) .and. (abs(gr%x(ip, idim)) <= bounds(1, idim)))
then
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
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, :) = &
1526 do ip_in = 1, bc%medium(idim)%points_number
1527 ip = bc%medium(idim)%points_map(ip_in)
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
1536 bc%medium(idim)%ep(ip_in) =
p_ep * (
m_one + ep_factor &
1538 bc%medium(idim)%mu(ip_in) =
p_mu * (
m_one + mu_factor &
1540 bc%medium(idim)%sigma_e(ip_in) = (
m_one + sigma_e_factor &
1542 bc%medium(idim)%sigma_m(ip_in) = (
m_one + sigma_m_factor &
1544 bc%medium(idim)%c(ip_in) =
m_one /
sqrt(bc%medium(idim)%ep(ip_in) * bc%medium(idim)%mu(ip_in))
1548 safe_deallocate_a(tmp)
1549 safe_deallocate_a(tmp_grad)
1558 class(
mesh_t),
intent(in) :: mesh
1560 real(real64),
intent(in) :: bounds(:,:)
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))
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))
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))
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))
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))
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))
1665 class(
mesh_t),
intent(in) :: mesh
1666 integer,
intent(in) :: ip
1667 real(real64),
intent(in) :: bounds(:,:)
1668 integer,
intent(out) :: point_info
1670 real(real64) :: rr, dd, xx(3), width(3)
1674 width(1:3) = bounds(2, 1:3) - bounds(1, 1:3)
1676 xx(1:mesh%box%dim) = mesh%x(ip, 1:mesh%box%dim)
1678 if (bc%ab_user_def)
then
1680 dd = bc%ab_ufn(ip) - bounds(1, 1)
1682 if (bc%ab_ufn(ip) < bounds(2, 1))
then
1689 select type (box => mesh%box)
1691 rr = norm2(xx - box%center)
1692 dd = rr - bounds(1, 1)
1694 if (dd < width(1))
then
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
1721 class(
mesh_t),
intent(in) :: mesh
1722 integer,
intent(in) :: ip
1723 real(real64),
intent(in) :: bounds(:,:)
1724 integer,
intent(out) :: boundary_info
1726 real(real64) :: xx(3)
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
1744 class(
mesh_t),
intent(in) :: mesh
1746 real(real64),
intent(in) :: bounds(:,:)
1748 integer :: ip, ip_in, ip_bd, point_info
1749 real(real64) :: xx(mesh%box%dim)
1760 if (all(abs(xx) <= bounds(2,1:mesh%box%dim)))
then
1761 if (any(abs(xx) > bounds(1,1:mesh%box%dim)))
then
1769 if (point_info == 0)
then
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.
1789 if (all(abs(xx) <= bounds(2,1:mesh%box%dim)))
then
1790 if (any(abs(xx) > bounds(1,1:mesh%box%dim)))
then
1798 if (point_info == 0)
then
1800 st%inner_points_map(ip_in) = ip
1801 st%inner_points_mask(ip) = .
true.
1804 st%boundary_points_map(ip_bd) = ip
1805 st%boundary_points_mask(ip) = .
true.
1811 call accel_write_buffer(st%buff_inner_points_map, st%inner_points_number, st%inner_points_map)
1813 call accel_write_buffer(st%buff_boundary_points_map, st%boundary_points_number, st%boundary_points_map)
1823 class(
mesh_t),
intent(in) :: mesh
1825 real(real64),
intent(in) :: bounds(:,:)
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)
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)
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)
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)
1851 safe_allocate(nn(1:2, 1:3, 1:3, 1:3))
1853 st%surface_grid_center(1, 1, :, :) = -int(bounds(1,1))
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))
1862 st%surface_grid_center(2, 1, :, :) = int(bounds(1,1))
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))
1872 st%surface_grid_center(1, 2, :, :) = -int(bounds(1,2))
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))
1881 st%surface_grid_center(2, 2, :, :) = int(bounds(1,2))
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))
1891 st%surface_grid_center(1, 3, :, :) = -int(bounds(1,3))
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))
1900 st%surface_grid_center(2, 3, :, :) = int(bounds(1,3))
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))
1910 st%surface_grid_points_number(:,:,:) = 0
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)
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)
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)
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)
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)
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)
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)
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)
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)
1981 safe_allocate(st%surface_grid_points_map(1:2, 1:st%dim, 1:ix_max, 1:iy_max, 1:nn_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)
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)
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))
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
2008 rr(2) = iiy * mesh%spacing(2)
2009 rr(3) = iiz * mesh%spacing(3)
2010 iix = int(bounds(1,1)/mesh%spacing(1))
2012 st%surface_grid_points_map(2, 1, idx1, idx2, nn(2, 1, idx1, idx2)) = ip_global
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)
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)
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))
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)
2041 rr(3) = iiz * mesh%spacing(3)
2042 iiy = int(bounds(1,2)/mesh%spacing(2))
2044 st%surface_grid_points_map(2, 2, idx1, idx2, nn(2, 2, idx1, idx2)) = ip_global
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)
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)
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))
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)
2074 iiz = int(bounds(1,3)/mesh%spacing(3))
2076 st%surface_grid_points_map(2, 3, idx1, idx2, nn(2, 3, idx1, idx2)) = ip_global
2081 safe_deallocate_a(nn)
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
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
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
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
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
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
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
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
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
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
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
subroutine, public accel_release_buffer(this, async)
pure logical function, public accel_is_enabled()
integer, parameter, public accel_mem_read_only
type(debug_t), save, public debug
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
real(real64), parameter, public m_huge
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public p_mu
real(real64), parameter, public p_ep
complex(real64), parameter, public m_z0
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
real(real64), parameter, public m_one
real(real64), parameter, public m_five
This module implements the underlying real-space grid.
This module implements the index, used for the mesh points.
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...
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.
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.
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.
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
character(len=512), private msg
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Some general things and nomenclature:
logical function, public parse_is_defined(namespace, name)
integer function, public parse_block(namespace, name, blk, check_varinfo_)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
type(type_t), public type_float
type(type_t), public type_cmplx
type(type_t), public type_integer
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
character(len=20) pure function, public units_abbrev(this)
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.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Describes mesh distribution to nodes.