66 integer,
parameter :: &
67 MEDIUM_PARALLELEPIPED = 1, &
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
81 type(multicomm_t) :: mc
82 type(single_medium_box_t) :: medium_box
98 procedure linear_medium_constructor
109 class(linear_medium_t),
pointer :: sys
110 type(namespace_t),
intent(in) :: namespace
128 class(linear_medium_t),
target,
intent(inout) :: this
129 type(namespace_t),
intent(in) :: namespace
131 integer :: nlines, ncols,n_points, n_global_points
132 integer(int64) :: index_range(4)
133 integer,
allocatable :: tmp(:)
140 this%namespace = namespace
141 this%space =
space_t(this%namespace)
142 if (this%space%is_periodic())
then
148 index_range(1) = this%gr%np_global
151 index_range(4) = 100000
155 index_range, (/ 5000, 1, 1, 1 /))
175 if (
parse_block(namespace,
'LinearMediumProperties', blk) == 0)
then
178 if (nlines /= 1)
then
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
199 message(1) =
'You must specify the properties of your linear medium through the LinearMediumProperties block.'
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'
223 allocate(this%supported_interactions(0))
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.))
235 this%outp%iter_dir, this%medium_box%points_map, this%medium_box%ep, &
236 this%medium_box%mu, this%medium_box%c)
238 if (this%medium_box%check_medium_points)
then
239 safe_allocate(tmp(1:this%gr%np))
241 write(
message(1),
'(a, a, a)')
'Check of points inside surface of medium ', trim(this%medium_box%filename),
":"
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
249 safe_deallocate_a(tmp)
264 select type (interaction)
266 message(1) =
"Trying to initialize an unsupported interaction by a linear medium."
283 select type (interaction)
286 partner%gr, partner%space, partner%namespace)
289 medium_box_grid => partner%medium_box%to_grid(partner%gr)
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)
300 safe_deallocate_p(regridding)
301 safe_deallocate_p(medium_box_grid)
303 message(1) =
"Trying to initialize an unsupported interaction by a linear medium."
323 character(len=:),
allocatable,
intent(out) :: updated_quantities(:)
336 real(real64),
intent(in) :: tol
354 select type (interaction)
358 message(1) =
"Unsupported interaction."
369 integer :: restart_file_unit
373 message(1) =
"Linear medium system "//trim(this%namespace%get())//
" will only write a dummy restart file."
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'
381 message(1) =
"Successfully wrote restart data for system "//trim(this%namespace%get())
392 integer :: restart_file_unit
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
407 message(1) =
"Successfully read restart data for system "//trim(this%namespace%get())
419 this%kinetic_energy =
m_zero
445 type(
grid_t),
intent(in) :: gr
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)
456 safe_allocate(tmp_points_map(1:gr%np))
457 safe_allocate(tmp_bdry_map(1:gr%np))
464 medium_box%global_points_number, tmp_points_map)
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
475 xx(1:3) = gr%x(ip, 1:3)
479 tmp_points_map(ip_in) = ip
483 tmp_bdry_map(ip_bd) = ip
487 medium_box%points_number = ip_in
488 medium_box%bdry_number = ip_bd
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)
496 medium_box%points_map(:) = tmp_points_map(1:medium_box%points_number)
499 safe_deallocate_a(tmp_points_map)
500 safe_deallocate_a(tmp_bdry_map)
507 logical pure function check_point_in_bounds(xx, bounds) result (check)
508 real(real64),
intent(in) :: xx(:)
509 real(real64),
intent(in) :: bounds(:,:)
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
520 logical pure function check_point_on_bounds(xx, bounds) result (check)
521 real(real64),
intent(in) :: xx(:)
522 real(real64),
intent(in) :: bounds(:,:)
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
548 type(single_medium_box_t),
intent(inout) :: medium_box
549 type(grid_t),
intent(in) :: gr
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(:,:)
557 call profiling_in(
'GET_LIN_MEDIUM_EM_PROP')
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))
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))
568 xx(1:3) = gr%x(ip,1:3)
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
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)))
588 else if (this%edge_profile == option__linearmediumedgeprofile__edged)
then
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
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)
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))
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)
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))
623 safe_deallocate_a(tmp)
624 safe_deallocate_a(tmp_grad)
626 call profiling_out(
'GET_LIN_MEDIUM_EM_PROP')
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
642 real(real64) :: xx(3)
643 type(cgal_polyhedra_t) :: cgal_poly
647 call profiling_in(
'GET_POINTS_MAP_FROM_FILE')
649 call cgal_polyhedron_init(cgal_poly, trim(filename), verbose = .false.)
653 if (
present(scale_factor))
then
654 xx(1:3) = scale_factor * mesh%x(ip, 1:3)
656 xx(1:3) = mesh%x(ip, 1:3)
658 if (cgal_polyhedron_point_inside(cgal_poly, xx(1), xx(2), xx(3)))
then
664 call cgal_polyhedron_end(cgal_poly)
666 call mpi_world%allreduce(ip_in, global_points_number, 1, mpi_integer, mpi_sum)
668 call profiling_out(
'GET_POINTS_MAP_FROM_FILE')
677 type(single_medium_box_t),
intent(inout) :: medium_box
678 type(namespace_t),
intent(in) :: namespace
680 integer :: nlines, ncols, idim
699 if (.not. varinfo_valid_option(
'LinearMediumBoxShape', medium_box%box_shape))
then
700 call messages_input_error(namespace,
'LinearMediumBoxShape')
703 select case (medium_box%box_shape)
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)
722 nlines = parse_block_n(blk)
723 if (nlines /= 1)
then
724 call messages_input_error(namespace,
'LinearMediumBoxSize',
'should consist of one line')
726 ncols = parse_block_cols(blk, 0)
728 call messages_input_error(namespace,
'LinearMediumBoxSize',
'should consist of six columns')
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))
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)
742 call messages_print_with_emphasis(namespace=namespace)
744 message(1) =
"For parallelepiped box shapes, you must provide a LinearMediumBoxSize block."
745 call messages_fatal(1, namespace=namespace)
755 if (parse_is_defined(namespace,
'LinearMediumBoxFile'))
then
756 call parse_variable(namespace,
'LinearMediumBoxFile',
'mediumboxfile', medium_box%filename)
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)
771 call parse_variable(namespace,
'CheckPointsMediumFromFile', .false., medium_box%check_medium_points)
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.
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
real(real64), parameter, public m_zero
character(len= *), parameter, public td_dir
This module implements the underlying real-space grid.
subroutine, public grid_init_stage_1(gr, namespace, space, symm, latt, n_sites, site_position)
First stage of the grid initialization.
subroutine, public grid_init_stage_2(gr, namespace, space, mc, qvector)
Second stage of the grid initialization.
subroutine, public grid_end(gr)
finalize a grid object
integer, parameter, public linear_medium_to_em_field
This module defines the abstract interaction_t class, and some auxiliary classes for interactions.
subroutine, public io_close(iunit, grp)
subroutine, public io_mkdir(fname, namespace, parents)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
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.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
subroutine, public messages_not_implemented(feature, namespace)
subroutine, public messages_warning(no_lines, all_nodes, namespace)
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)
type(mpi_grp_t), public mpi_world
This module handles the communicators for the various parallelization strategies.
subroutine, public multicomm_end(mc)
subroutine, public multicomm_init(mc, namespace, base_grp, mode_para, n_node, index_range, min_range)
create index and domain communicators
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
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.
This module implements the basic propagator framework.
This module defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
Implementation details for regridding.
This module implements the abstract system type.
subroutine, public system_end(this)
Descriptor of one algorithmic operation.
Description of the grid, containing information on derivatives, stencil, and symmetries.
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...
contains the information of the meshes and provides the transfer functions
Abstract class for systems.