Octopus
linear_medium.F90
Go to the documentation of this file.
1!! Copyright (C) 2020 F. Bonafé, H. Appel, R. Jestädt
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
19#include "global.h"
20
27 use comm_oct_m
28 use debug_oct_m
30 use global_oct_m
31 use grid_oct_m
35 use iso_c_binding
36 use io_oct_m
37 use math_oct_m
39 use math_oct_m
41 use mesh_oct_m
42 use mpi_oct_m
47 use parser_oct_m
53 use space_oct_m
54 use system_oct_m
56
57 implicit none
58
59 private
60 public :: &
65
66 integer, parameter :: &
67 MEDIUM_PARALLELEPIPED = 1, &
69
72 type, extends(system_t) :: linear_medium_t
73 type(space_t) :: space
74 real(real64) :: ep_factor
75 real(real64) :: mu_factor
76 real(real64) :: sigma_e_factor
77 real(real64) :: sigma_m_factor
78 integer :: edge_profile
79 type(output_t) :: outp
80 type(grid_t) :: gr
81 type(multicomm_t) :: mc
82 type(single_medium_box_t) :: medium_box
83
84 contains
85 procedure :: init_interaction => linear_medium_init_interaction
86 procedure :: init_interaction_as_partner => linear_medium_init_interaction_as_partner
87 procedure :: initialize => linear_medium_initialize
88 procedure :: do_algorithmic_operation => linear_medium_do_algorithmic_operation
89 procedure :: is_tolerance_reached => linear_medium_is_tolerance_reached
90 procedure :: copy_quantities_to_interaction => linear_medium_copy_quantities_to_interaction
91 procedure :: restart_write_data => linear_medium_restart_write_data
92 procedure :: restart_read_data => linear_medium_restart_read_data
93 procedure :: update_kinetic_energy => linear_medium_update_kinetic_energy
95 end type linear_medium_t
96
97 interface linear_medium_t
98 procedure linear_medium_constructor
99 end interface linear_medium_t
100
101contains
102
103 ! ---------------------------------------------------------
108 function linear_medium_constructor(namespace) result(sys)
109 class(linear_medium_t), pointer :: sys
110 type(namespace_t), intent(in) :: namespace
111
113
114 allocate(sys)
115
116 call linear_medium_init(sys, namespace)
117
119 end function linear_medium_constructor
120
121 ! ---------------------------------------------------------
126 ! ---------------------------------------------------------
127 subroutine linear_medium_init(this, namespace)
128 class(linear_medium_t), target, intent(inout) :: this
129 type(namespace_t), intent(in) :: namespace
130
131 integer :: nlines, ncols,n_points, n_global_points
132 integer(int64) :: index_range(4)
133 integer, allocatable :: tmp(:)
134 type(block_t) :: blk
135
136 push_sub(linear_medium_init)
137
138 call profiling_in('MEDIUM_BOX_INIT')
139
140 this%namespace = namespace
141 this%space = space_t(this%namespace)
142 if (this%space%is_periodic()) then
143 call messages_not_implemented('Linear medium for periodic systems', namespace=namespace)
144 end if
145 call grid_init_stage_1(this%gr, this%namespace, this%space)
146 ! store the ranges for these two indices (serves as initial guess
147 ! for parallelization strategy)
148 index_range(1) = this%gr%np_global ! Number of points in mesh
149 index_range(2) = 1 ! Number of states
150 index_range(3) = 1 ! Number of k-points
151 index_range(4) = 100000 ! Some large number
152
153 ! create index and domain communicators
154 call multicomm_init(this%mc, this%namespace, mpi_world, calc_mode_par, mpi_world%size, &
155 index_range, (/ 5000, 1, 1, 1 /))
156 call grid_init_stage_2(this%gr, this%namespace, this%space, this%mc)
157
158 call medium_box_init(this%medium_box, this%namespace)
160 !%Variable LinearMediumProperties
161 !%Type block
162 !%Section Maxwell::Medium
163 !%Description
164 !% Defines electromagnetic parameters for a linear medium box.
165 !%
166 !% Example:
167 !%
168 !% <tt>%LinearMediumProperties
169 !% <br>&nbsp;&nbsp; epsilon_factor | mu_factor | sigma_e | sigma_m
170 !% <br>%</tt>
171 !%
172 !% Permittivity factor, permeability factor, electric conductivity and magnetic conductivity of the medium box.
173 !%End
175 if (parse_block(namespace, 'LinearMediumProperties', blk) == 0) then
176
177 nlines = parse_block_n(blk)
178 if (nlines /= 1) then
179 call messages_input_error(namespace, 'LinearMediumProperties', 'should consist of one line', nlines)
180 end if
181 ncols = parse_block_cols(blk, 0)
182 if (ncols /= 4) then
183 call messages_input_error(namespace, 'LinearMediumProperties', 'should consist of four columns')
184 end if
185 call parse_block_float(blk, 0, 0, this%ep_factor)
186 call parse_block_float(blk, 0, 1, this%mu_factor)
187 call parse_block_float(blk, 0, 2, this%sigma_e_factor)
188 call parse_block_float(blk, 0, 3, this%sigma_m_factor)
189 write(message(1),'(a,es9.2)') 'Box epsilon factor: ', this%ep_factor
190 write(message(2),'(a,es9.2)') 'Box mu factor: ', this%mu_factor
191 write(message(3),'(a,es9.2)') 'Box electric sigma: ', this%sigma_e_factor
192 write(message(4),'(a,es9.2)') 'Box magnetic sigma: ', this%sigma_m_factor
193 write(message(5),'(a)') ""
194 call messages_info(5, namespace=namespace)
195 call parse_block_end(blk)
196
197 call messages_print_with_emphasis(namespace=namespace)
198 else
199 message(1) = 'You must specify the properties of your linear medium through the LinearMediumProperties block.'
200 call messages_fatal(1, namespace=namespace)
201 end if
202
203 !%Variable LinearMediumEdgeProfile
204 !%Type integer
205 !%Section Maxwell::Medium
206 !%Description
207 !% Defines the type of numerical approximation used for the derivatives at the edges of the medium box.
208 !% Default is edged. When the box shape is read from file, only the edged profile is supported.
209 !%
210 !%Option edged 1
211 !% Medium box edges are considered steep for derivatives.
212 !%Option smooth 2
213 !% Medium box edged and softened for derivatives.
214 !%End
215 call parse_variable(namespace, 'LinearMediumEdgeProfile', option__linearmediumedgeprofile__edged, this%edge_profile)
216 if (this%edge_profile == option__linearmediumedgeprofile__edged) then
217 write(message(1),'(a,a)') 'Box shape: ', 'edged'
218 else if (this%edge_profile == option__linearmediumedgeprofile__smooth) then
219 write(message(1),'(a,a)') 'Box shape: ', 'smooth'
220 end if
221 call messages_info(1, namespace=namespace)
222
223 allocate(this%supported_interactions(0))
224 this%supported_interactions_as_partner = [linear_medium_to_em_field]
225 call this%quantities%add(quantity_t("permitivity", updated_on_demand = .false., always_available = .true.))
226 call this%quantities%add(quantity_t("permeability", updated_on_demand = .false., always_available = .true.))
227 call this%quantities%add(quantity_t("E conductivity", updated_on_demand = .false., always_available = .true.))
228 call this%quantities%add(quantity_t("M conductivity", updated_on_demand = .false., always_available = .true.))
229
230 call output_linear_medium_init(this%outp, this%namespace, this%space)
231
232 call get_medium_box_points_map(this%medium_box, this%gr)
233 call get_linear_medium_em_properties(this, this%medium_box, this%gr)
234 call output_linear_medium(this%outp, this%namespace, this%space, this%gr, &
235 this%outp%iter_dir, this%medium_box%points_map, this%medium_box%ep, &
236 this%medium_box%mu, this%medium_box%c)
237
238 if (this%medium_box%check_medium_points) then
239 safe_allocate(tmp(1:this%gr%np))
240 n_global_points = 0
241 write(message(1),'(a, a, a)') 'Check of points inside surface of medium ', trim(this%medium_box%filename), ":"
242 call messages_info(1, namespace=this%namespace)
243 call get_points_map_from_file(this%medium_box%filename, this%gr, n_points, n_global_points, tmp, 0.99_real64)
244 write(message(1),'(a, I8)')'Number of points inside medium (normal coordinates):', &
245 this%medium_box%global_points_number
246 write(message(2),'(a, I8)')'Number of points inside medium (rescaled coordinates):', n_global_points
247 write(message(3), '(a)') ""
248 call messages_info(3, namespace=this%namespace)
249 safe_deallocate_a(tmp)
250 end if
251
252 call profiling_out('MEDIUM_BOX_INIT')
253
254 pop_sub(linear_medium_init)
255 end subroutine linear_medium_init
256
257 ! ---------------------------------------------------------
258 subroutine linear_medium_init_interaction(this, interaction)
259 class(linear_medium_t), target, intent(inout) :: this
260 class(interaction_t), intent(inout) :: interaction
261
263
264 select type (interaction)
265 class default
266 message(1) = "Trying to initialize an unsupported interaction by a linear medium."
267 call messages_fatal(1, namespace=this%namespace)
268 end select
269
271 end subroutine linear_medium_init_interaction
272
273 ! ---------------------------------------------------------
274 subroutine linear_medium_init_interaction_as_partner(partner, interaction)
275 class(linear_medium_t), intent(in) :: partner
276 class(interaction_surrogate_t), intent(inout) :: interaction
277
278 type(regridding_t), pointer :: regridding
279 type(single_medium_box_t), pointer :: medium_box_grid
280
282
283 select type (interaction)
285 regridding => regridding_t(interaction%system_gr, &
286 partner%gr, partner%space, partner%namespace)
287
288 ! create a medium box with the size of the partner grid and map all quantities to it
289 medium_box_grid => partner%medium_box%to_grid(partner%gr)
290
291 ! now do the grid transfer to the medium box stored in the interaction which has the size of the system grid
292 call regridding%do_transfer(interaction%medium_box%aux_ep, medium_box_grid%aux_ep)
293 call regridding%do_transfer(interaction%medium_box%aux_mu, medium_box_grid%aux_mu)
294 call regridding%do_transfer(interaction%medium_box%ep, medium_box_grid%ep)
295 call regridding%do_transfer(interaction%medium_box%mu, medium_box_grid%mu)
296 call regridding%do_transfer(interaction%medium_box%c, medium_box_grid%c)
297 call regridding%do_transfer(interaction%medium_box%sigma_e, medium_box_grid%sigma_e)
298 call regridding%do_transfer(interaction%medium_box%sigma_m, medium_box_grid%sigma_m)
299
300 safe_deallocate_p(regridding)
301 safe_deallocate_p(medium_box_grid)
302 class default
303 message(1) = "Trying to initialize an unsupported interaction by a linear medium."
304 call messages_fatal(1, namespace=partner%namespace)
305 end select
306
309
310 ! ---------------------------------------------------------
311 subroutine linear_medium_initialize(this)
312 class(linear_medium_t), intent(inout) :: this
313
315
317 end subroutine linear_medium_initialize
318
319 ! ---------------------------------------------------------
320 logical function linear_medium_do_algorithmic_operation(this, operation, updated_quantities) result(done)
321 class(linear_medium_t), intent(inout) :: this
322 class(algorithmic_operation_t), intent(in) :: operation
323 character(len=:), allocatable, intent(out) :: updated_quantities(:)
324
326
327 ! Currently there are no linear medium specific algorithmic operations
328 done = .false.
329
332
333 ! ---------------------------------------------------------
334 logical function linear_medium_is_tolerance_reached(this, tol) result(converged)
335 class(linear_medium_t), intent(in) :: this
336 real(real64), intent(in) :: tol
337
339
340 ! this routine is never called at present, no reason to be here
341 assert(.false.)
342 converged = .false.
343
346
347 ! ---------------------------------------------------------
348 subroutine linear_medium_copy_quantities_to_interaction(partner, interaction)
349 class(linear_medium_t), intent(inout) :: partner
350 class(interaction_surrogate_t), intent(inout) :: interaction
353
354 select type (interaction)
356 ! No need to copy anything
357 class default
358 message(1) = "Unsupported interaction."
359 call messages_fatal(1, namespace=partner%namespace)
360 end select
361
364
365 ! ---------------------------------------------------------
366 subroutine linear_medium_restart_write_data(this)
367 class(linear_medium_t), intent(inout) :: this
368
369 integer :: restart_file_unit
370
372
373 message(1) = "Linear medium system "//trim(this%namespace%get())//" will only write a dummy restart file."
374 call messages_warning(1, namespace=this%namespace)
375
376 call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
377 restart_file_unit = io_open('restart/'//td_dir// 'restart_linear_medium', this%namespace, action='write')
378 write(restart_file_unit, *) 'Linear medium restart file'
379 call io_close(restart_file_unit)
380
381 message(1) = "Successfully wrote restart data for system "//trim(this%namespace%get())
382 call messages_info(1, namespace=this%namespace)
383
386
387 ! ---------------------------------------------------------
388 ! this function returns true if restart data could be read
389 logical function linear_medium_restart_read_data(this)
390 class(linear_medium_t), intent(inout) :: this
391
392 integer :: restart_file_unit
393
395
396 call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
397 restart_file_unit = io_open('restart/'//td_dir// 'restart_linear_medium', this%namespace, action='read', die=.false.)
398 if (restart_file_unit /= -1) then
399 call io_close(restart_file_unit)
401 else
402 ! could not open file
404 end if
405
407 message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
408 call messages_info(1, namespace=this%namespace)
409 end if
410
415 class(linear_medium_t), intent(inout) :: this
416
418
419 this%kinetic_energy = m_zero
420
422
424
425 ! ---------------------------------------------------------
426 subroutine linear_medium_finalize(this)
427 type(linear_medium_t), intent(inout) :: this
428
429 push_sub(linear_medium_finalize)
430
431 call single_medium_box_end(this%medium_box)
432 call system_end(this)
433 call multicomm_end(this%mc)
434 call grid_end(this%gr)
435
437 end subroutine linear_medium_finalize
438
439
440 ! Specific routines for this system:
442 ! ---------------------------------------------------------
443 subroutine get_medium_box_points_map(medium_box, gr)
444 type(single_medium_box_t), intent(inout) :: medium_box
445 type(grid_t), intent(in) :: gr
446
447 integer :: ip, ip_in, ip_bd, idim
448 integer, allocatable :: tmp_points_map(:), tmp_bdry_map(:)
449 real(real64) :: bounds(2,gr%box%dim), xx(gr%box%dim)
450 logical :: inside
451
453
454 call profiling_in('GET_MEDIUM_BOX_POINTS_MAP')
455
456 safe_allocate(tmp_points_map(1:gr%np))
457 safe_allocate(tmp_bdry_map(1:gr%np))
458 tmp_points_map = 0
459 tmp_bdry_map = 0
460
461 if (medium_box%box_shape == medium_box_file) then
462
463 call get_points_map_from_file(medium_box%filename, gr, medium_box%points_number,&
464 medium_box%global_points_number, tmp_points_map)
465
466 else
467
468 do idim = 1, 3
469 bounds(1,idim) = medium_box%center(idim) - medium_box%lsize(idim)/m_two
470 bounds(2,idim) = medium_box%center(idim) + medium_box%lsize(idim)/m_two
471 end do
472 ip_in = 0
473 ip_bd = 0
474 do ip = 1, gr%np
475 xx(1:3) = gr%x(ip, 1:3)
476 inside = check_point_in_bounds(xx, bounds(:,:))
477 if (check_point_in_bounds(xx, bounds(:,:))) then
478 ip_in = ip_in + 1
479 tmp_points_map(ip_in) = ip
480 end if
481 if (check_point_on_bounds(xx, bounds(:,:))) then
482 ip_bd = ip_bd + 1
483 tmp_bdry_map(ip_bd) = ip
484 end if
485 end do
486
487 medium_box%points_number = ip_in
488 medium_box%bdry_number = ip_bd
489
490 safe_allocate(medium_box%bdry_map(1:medium_box%bdry_number))
491 medium_box%bdry_map = 0
492 medium_box%bdry_map = tmp_bdry_map(1:medium_box%bdry_number)
493
494 end if
495 call single_medium_box_allocate(medium_box, medium_box%points_number)
496 medium_box%points_map(:) = tmp_points_map(1:medium_box%points_number)
497
498 ! TODO: add some check that medium boxes do not overlap (now they are different systems)
499 safe_deallocate_a(tmp_points_map)
500 safe_deallocate_a(tmp_bdry_map)
501
502 call profiling_out('GET_MEDIUM_BOX_POINTS_MAP')
503
505 contains
506
507 logical pure function check_point_in_bounds(xx, bounds) result (check)
508 real(real64), intent(in) :: xx(:)
509 real(real64), intent(in) :: bounds(:,:)
510
511 check = .false.
512 if ((xx(1) >= bounds(1,1)) .and. (xx(1) <= bounds(2,1)) .and. &
513 (xx(2) >= bounds(1,2)) .and. (xx(2) <= bounds(2,2)) .and. &
514 (xx(3) >= bounds(1,3)) .and. (xx(3) <= bounds(2,3)) ) then
515 check = .true.
516 end if
517
518 end function check_point_in_bounds
520 logical pure function check_point_on_bounds(xx, bounds) result (check)
521 real(real64), intent(in) :: xx(:)
522 real(real64), intent(in) :: bounds(:,:)
523
524 check = .false.
525 if (is_close(xx(1), bounds(1,1)) .and. (xx(2) >= bounds(1,2) .and. xx(3) >= bounds(1,3)) &
526 .and. (xx(2) <= bounds(2,2) .and. xx(3) <= bounds(2,3)) .or. &
527 is_close(xx(2), bounds(1,2)) .and. (xx(1) >= bounds(1,1) .and. xx(3) >= bounds(1,3)) &
528 .and. (xx(1) <= bounds(2,1) .and. xx(3) <= bounds(2,3)) .or. &
529 is_close(xx(3), bounds(1,3)) .and. (xx(1) >= bounds(1,1) .and. xx(2) >= bounds(1,2)) &
530 .and. (xx(1) <= bounds(2,1) .and. xx(2) <= bounds(2,2)) .or. &
531 is_close(xx(1), bounds(2,1)) .and. (xx(2) >= bounds(1,2) .and. xx(3) >= bounds(1,3)) &
532 .and. (xx(2) <= bounds(2,2) .and. xx(3) <= bounds(2,3)) .or. &
533 is_close(xx(2), bounds(2,2)) .and. (xx(1) >= bounds(1,1) .and. xx(3) >= bounds(1,3)) &
534 .and. (xx(1) <= bounds(2,1) .and. xx(3) <= bounds(2,3)) .or. &
535 is_close(xx(3), bounds(2,3)) .and. (xx(1) >= bounds(1,1) .and. xx(2) >= bounds(1,2)) &
536 .and. (xx(1) <= bounds(2,1) .and. xx(2) <= bounds(2,2)) ) then
537 check = .true.
538 end if
539
540 end function check_point_on_bounds
541
542 end subroutine get_medium_box_points_map
543
544 ! ----------------------------------------------------------
546 subroutine get_linear_medium_em_properties(this, medium_box, gr)
547 type(linear_medium_t), intent(in) :: this
548 type(single_medium_box_t), intent(inout) :: medium_box
549 type(grid_t), intent(in) :: gr
550
551 integer :: ip, ip_in, ip_bd, ipp
552 real(real64) :: xx(gr%box%dim), xxp(gr%box%dim), dd, dd_max, dd_min
553 real(real64), allocatable :: tmp(:), tmp_grad(:,:)
554
556
557 call profiling_in('GET_LIN_MEDIUM_EM_PROP')
558
559 safe_allocate(tmp(1:gr%np_part))
560 safe_allocate(tmp_grad(1:gr%np_part, 1:gr%box%dim))
561 dd_max = max(2*gr%spacing(1), 2*gr%spacing(2), 2*gr%spacing(3))
562
563 do ip_in = 1, medium_box%points_number
564 ip = medium_box%points_map(ip_in)
565 if (this%edge_profile == option__linearmediumedgeprofile__smooth) then
566 assert(allocated(medium_box%bdry_map))
567
568 xx(1:3) = gr%x(ip,1:3)
569 dd_min = m_huge
570
571 do ip_bd = 1, medium_box%bdry_number
572 ipp = medium_box%bdry_map(ip_bd)
573 xxp(1:3) = gr%x(ipp,1:3)
574 dd = norm2(xx(1:3) - xxp(1:3))
575 if (dd < dd_min) dd_min = dd
576 end do
577
578 medium_box%ep(ip_in) = p_ep + ((p_ep * this%ep_factor - p_ep) &
579 * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max))))
580 medium_box%mu(ip_in) = p_mu + ((p_mu * this%mu_factor - p_mu) &
581 * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max))))
582 medium_box%c(ip_in) = m_one/sqrt(medium_box%ep(ip_in)*medium_box%mu(ip_in))
583 medium_box%sigma_e(ip_in) = this%sigma_e_factor &
584 * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max)))
585 medium_box%sigma_m(ip_in) = this%sigma_m_factor &
586 * m_one/(m_one + exp(-m_five/dd_max * (dd_min - m_two*dd_max)))
587
588 else if (this%edge_profile == option__linearmediumedgeprofile__edged) then
589
590 medium_box%ep(ip_in) = p_ep * this%ep_factor
591 medium_box%mu(ip_in) = p_mu * this%mu_factor
592 medium_box%c(ip_in) = m_one/sqrt(medium_box%ep(ip_in)*medium_box%mu(ip_in))
593 medium_box%sigma_e(ip_in) = this%sigma_e_factor
594 medium_box%sigma_m(ip_in) = this%sigma_m_factor
595
596 end if
597 end do
598
599 tmp(:) = p_ep
600 do ip_in = 1, medium_box%points_number
601 ip = medium_box%points_map(ip_in)
602 tmp(ip)= medium_box%ep(ip_in)
603 end do
604 call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
605 do ip_in = 1, medium_box%points_number
606 ip = medium_box%points_map(ip_in)
607 medium_box%aux_ep(ip_in, :) = tmp_grad(ip, :)/(m_four * medium_box%ep(ip_in))
608 end do
609
610 tmp(:) = p_mu
611 do ip_in = 1, medium_box%points_number
612 ip = medium_box%points_map(ip_in)
613 tmp(ip) = medium_box%mu(ip_in)
614 end do
615 call dderivatives_grad(gr%der, tmp, tmp_grad, set_bc = .false.)
616 do ip_in = 1, medium_box%points_number
617 ip = medium_box%points_map(ip_in)
618 medium_box%aux_mu(ip_in, :) = tmp_grad(ip, :)/(m_four * medium_box%mu(ip_in))
619 end do
620
621 !TODO: add print information about the medium box
622
623 safe_deallocate_a(tmp)
624 safe_deallocate_a(tmp_grad)
625
626 call profiling_out('GET_LIN_MEDIUM_EM_PROP')
627
630
631 ! ----------------------------------------------------------
633 subroutine get_points_map_from_file(filename, mesh, n_points, global_points_number, tmp_map, scale_factor)
634 character(len=256), intent(in) :: filename
635 class(mesh_t), intent(in) :: mesh
636 integer, intent(out) :: n_points
637 integer, intent(out) :: global_points_number
638 integer, intent(inout) :: tmp_map(:)
639 real(real64), optional, intent(in) :: scale_factor
640
641 integer :: ip_in, ip
642 real(real64) :: xx(3)
643 type(cgal_polyhedra_t) :: cgal_poly
644
646
647 call profiling_in('GET_POINTS_MAP_FROM_FILE')
648
649 call cgal_polyhedron_init(cgal_poly, trim(filename), verbose = .false.)
650
651 ip_in = 0
652 do ip = 1, mesh%np
653 if (present(scale_factor)) then
654 xx(1:3) = scale_factor * mesh%x(ip, 1:3)
655 else
656 xx(1:3) = mesh%x(ip, 1:3)
657 end if
658 if (cgal_polyhedron_point_inside(cgal_poly, xx(1), xx(2), xx(3))) then
659 ip_in = ip_in + 1
660 tmp_map(ip_in) = ip
661 end if
662 end do
663 n_points = ip_in
664 call cgal_polyhedron_end(cgal_poly)
665
666 call mpi_world%allreduce(ip_in, global_points_number, 1, mpi_integer, mpi_sum)
667
668 call profiling_out('GET_POINTS_MAP_FROM_FILE')
669
671
672 end subroutine get_points_map_from_file
673
674 ! ----------------------------------------------------------
676 subroutine medium_box_init(medium_box, namespace)
677 type(single_medium_box_t), intent(inout) :: medium_box
678 type(namespace_t), intent(in) :: namespace
679
680 integer :: nlines, ncols, idim
681 type(block_t) :: blk
682
683 push_sub(medium_box_init)
684
685 !%Variable LinearMediumBoxShape
686 !%Type integer
687 !%Section Maxwell::Medium
688 !%Description
689 !% This variable defines the shape of the linear medium box.
690 !% The default is <tt>medium_parallelepiped</tt>.
691 !%Option medium_parallelepiped 1
692 !% The medium box will be a parallelepiped whose center and dimensions are taken from
693 !% the variable <tt>LinearMediumBoxSize</tt>.
694 !%Option medium_box_file 2
695 !% The simulation box will be read from an external file in OFF format, defined by the variable <tt>LinearMediumBoxFile</tt>.
696 !%End
697 call parse_variable(namespace, 'LinearMediumBoxShape', medium_parallelepiped, medium_box%box_shape)
698
699 if (.not. varinfo_valid_option('LinearMediumBoxShape', medium_box%box_shape)) then
700 call messages_input_error(namespace, 'LinearMediumBoxShape')
701 end if
702
703 select case (medium_box%box_shape)
705
706 !%Variable LinearMediumBoxSize
707 !%Type block
708 !%Section Maxwell::Medium
709 !%Description
710 !% Defines center and size of a parallelepiped linear medium box.
711 !%
712 !% Example:
713 !%
714 !% <tt>%LinearMediumBoxSize
715 !% <br>&nbsp;&nbsp; center_x | center_y | center_z | x_length | y_length | z_length
716 !% <br>%</tt>
717 !%End
718
719 if (parse_block(namespace, 'LinearMediumBoxSize', blk) == 0) then
720 call messages_print_with_emphasis(msg=trim('Linear Medium box center and size:'), namespace=namespace)
721
722 nlines = parse_block_n(blk)
723 if (nlines /= 1) then
724 call messages_input_error(namespace, 'LinearMediumBoxSize', 'should consist of one line')
725 end if
726 ncols = parse_block_cols(blk, 0)
727 if (ncols /= 6) then
728 call messages_input_error(namespace, 'LinearMediumBoxSize', 'should consist of six columns')
729 end if
730 do idim = 1, 3
731 call parse_block_float(blk, 0, idim-1, medium_box%center(idim))
732 call parse_block_float(blk, 0, idim+2, medium_box%lsize(idim))
733 end do
734 write(message(1),'(a,es9.2,a,es9.2,a,es9.2)') 'Box center: ', medium_box%center(1), ' | ',&
735 medium_box%center(2), ' | ', medium_box%center(3)
736 write(message(2),'(a,es9.2,a,es9.2,a,es9.2)') 'Box size : ', medium_box%lsize(1), ' | ', &
737 medium_box%lsize(2), ' | ', medium_box%lsize(3)
738 write(message(3),'(a)') ""
739 call messages_info(3)
740 call parse_block_end(blk)
741
742 call messages_print_with_emphasis(namespace=namespace)
743 else
744 message(1) = "For parallelepiped box shapes, you must provide a LinearMediumBoxSize block."
745 call messages_fatal(1, namespace=namespace)
746 end if
747
748 case (medium_box_file)
749 !%Variable LinearMediumBoxFile
750 !%Type string
751 !%Section Maxwell::Medium
752 !%Description
753 !% File in OFF format with the shape of the linear medium.
754 !%End
755 if (parse_is_defined(namespace, 'LinearMediumBoxFile')) then
756 call parse_variable(namespace, 'LinearMediumBoxFile', 'mediumboxfile', medium_box%filename)
757 else
758 message(1) = "When using box_file as the box shape, you must provide a filename through the LinearMediumBoxFile variable."
759 call messages_fatal(1, namespace=namespace)
760 end if
761
762 !%Variable CheckPointsMediumFromFile
763 !%Type logical
764 !%Default no
765 !%Section Maxwell::Medium
766 !%Description
767 !% Whether to re-calculate the points map by artificially shrinking the coordinate system by a factor of
768 !% 0.99 to check if the points inside the medium surface are properly detected. This works for only one
769 !% medium surface which is centered in the origin of the coordinate system.
770 !%End
771 call parse_variable(namespace, 'CheckPointsMediumFromFile', .false., medium_box%check_medium_points)
772
773 end select
774
775 pop_sub(medium_box_init)
776
777 end subroutine medium_box_init
778
779end module linear_medium_oct_m
780
781!! Local Variables:
782!! mode: f90
783!! coding: utf-8
784!! End:
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
logical pure function check_point_on_bounds(xx, bounds)
logical pure function check_point_in_bounds(xx, bounds)
This module implements the basic elements defining algorithms.
Definition: algorithm.F90:141
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
character(len= *), parameter, public td_dir
Definition: global.F90:251
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public grid_init_stage_1(gr, namespace, space, symm, latt, n_sites, site_position)
First stage of the grid initialization.
Definition: grid.F90:194
subroutine, public grid_init_stage_2(gr, namespace, space, mc, qvector)
Second stage of the grid initialization.
Definition: grid.F90:465
subroutine, public grid_end(gr)
finalize a grid object
Definition: grid.F90:493
integer, parameter, public linear_medium_to_em_field
This module defines the abstract interaction_t class, and some auxiliary classes for interactions.
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:311
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
This module defines a linear medium for use in classical electrodynamics calculations.
subroutine linear_medium_initialize(this)
logical function linear_medium_is_tolerance_reached(this, tol)
integer, parameter medium_parallelepiped
subroutine linear_medium_copy_quantities_to_interaction(partner, interaction)
subroutine linear_medium_restart_write_data(this)
class(linear_medium_t) function, pointer linear_medium_constructor(namespace)
The factory routine (or constructor) allocates a pointer of the corresponding type and then calls the...
logical function linear_medium_restart_read_data(this)
subroutine, public get_medium_box_points_map(medium_box, gr)
subroutine linear_medium_init_interaction(this, interaction)
subroutine, public medium_box_init(medium_box, namespace)
Parse and store geometry of medium box.
subroutine linear_medium_init_interaction_as_partner(partner, interaction)
subroutine linear_medium_update_kinetic_energy(this)
integer, parameter medium_box_file
logical function linear_medium_do_algorithmic_operation(this, operation, updated_quantities)
subroutine linear_medium_finalize(this)
subroutine get_points_map_from_file(filename, mesh, n_points, global_points_number, tmp_map, scale_factor)
Populate list of point indices for points inside the polyhedron.
subroutine, public linear_medium_init(this, namespace)
The init routine is a module level procedure This has the advantage that different classes can have d...
subroutine get_linear_medium_em_properties(this, medium_box, gr)
Evaluate electromagnetic properties of linear medium.
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:115
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
subroutine, public multicomm_end(mc)
Definition: multicomm.F90:889
subroutine, public multicomm_init(mc, namespace, base_grp, mode_para, n_node, index_range, min_range)
create index and domain communicators
Definition: multicomm.F90:264
subroutine, public output_linear_medium(outp, namespace, space, mesh, dir, points_map, ep, mu, cc)
subroutine, public output_linear_medium_init(outp, namespace, space)
this module contains the low-level part of the output system
Definition: output_low.F90:115
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
This module implements the basic propagator framework.
Definition: propagator.F90:117
This module defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
Definition: quantity.F90:138
Implementation details for regridding.
Definition: regridding.F90:170
This module implements the abstract system type.
Definition: system.F90:118
subroutine, public system_end(this)
Definition: system.F90:1145
Descriptor of one algorithmic operation.
Definition: algorithm.F90:163
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:168
abstract interaction class
surrogate interaction class to avoid circular dependencies between modules.
linear medium for classical electrodynamics
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...
Definition: quantity.F90:171
contains the information of the meshes and provides the transfer functions
Definition: regridding.F90:197
Abstract class for systems.
Definition: system.F90:172
int true(void)