37 use,
intrinsic :: iso_fortran_env
72 type(space_t) :: space
73 real(real64) :: omega_p
74 real(real64) :: gamma_p
75 real(real64) :: strength_p
76 real(real64),
allocatable :: current_p(:,:)
77 real(real64),
allocatable :: e_field(:,:)
78 real(real64),
allocatable :: e_field_dt_half(:,:)
79 real(real64),
allocatable :: e_field_dt_full(:,:)
80 real(real64),
allocatable :: current_at_point(:,:)
81 real(real64),
allocatable :: selected_points_coordinate(:,:)
82 integer :: n_output_points
83 integer :: medium_type
85 type(multicomm_t) :: mc
86 type(c_ptr) :: write_handle
87 type(output_t) :: outp
88 logical :: from_scratch = .
true.
89 type(restart_t) :: restart_load
90 type(restart_t) :: restart_dump
111 procedure dispersive_medium_constructor
114 integer,
public,
parameter :: &
125 class(dispersive_medium_t),
pointer :: sys
126 type(namespace_t),
intent(in) :: namespace
144 class(dispersive_medium_t),
target,
intent(inout) :: this
145 type(namespace_t),
intent(in) :: namespace
147 integer :: nlines, ncols, idim, il
148 real(real64) :: pos(3)
153 this%namespace = namespace
157 this%space =
space_t(this%namespace)
158 if (this%space%is_periodic())
then
161 if (this%space%dim /= 3)
then
177 call parse_variable(namespace,
'MediumDispersionType', drude_medium, this%medium_type)
223 if (
parse_block(namespace,
'MediumCurrentCoordinates', blk) == 0)
then
225 this%n_output_points = nlines
226 safe_allocate(this%selected_points_coordinate(1:nlines,1:3))
227 safe_allocate(this%current_at_point(1:nlines,1:3))
233 this%selected_points_coordinate(il,:) = pos(:)
234 this%current_at_point(il,:) =
m_zero
238 this%n_output_points = 1
239 safe_allocate(this%selected_points_coordinate(1,3))
240 safe_allocate(this%current_at_point(1,3))
241 this%selected_points_coordinate =
m_zero
242 this%current_at_point =
m_zero
248 call this%quantities%add(
quantity_t(
"current", updated_on_demand = .false.))
260 integer(int64) :: index_range(4)
268 index_range(1) = this%gr%np_global
271 index_range(4) = 100000
275 index_range, (/ 5000, 1, 1, 1 /))
279 this%mc, ierr, mesh=this%gr)
281 this%mc, ierr, mesh=this%gr)
283 safe_allocate(this%current_p(1:this%gr%np, 1:3))
284 safe_allocate(this%e_field(1:this%gr%np, 1:3))
285 this%e_field(:,:) =
m_zero
299 select type (interaction)
301 call interaction%init(this%gr, 3)
306 select type (prop => this%algo)
311 message(1) =
"The chosen propagator does not yet support interaction interpolation"
314 call interaction%init_interpolation(depth, interaction%label)
317 message(1) =
"Trying to initialize an unsupported interaction by a Dispersive medium."
331 select type (interaction)
333 call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
335 message(1) =
"Trying to initialize an unsupported interaction by a linear medium."
348 this%from_scratch = .
true.
358 character(len=:),
allocatable,
intent(out) :: updated_quantities(:)
361 real(real64) :: k1(1:3), k2(1:3), k3(1:3), k4(1:3)
364 call profiling_in(trim(this%namespace%get())//
":"//trim(operation%id))
374 select type (algo => this%algo)
378 select case (operation%id)
383 safe_allocate(this%e_field_dt_half(1:this%gr%np, 1:3))
384 safe_allocate(this%e_field_dt_full(1:this%gr%np, 1:3))
387 safe_deallocate_a(this%e_field_dt_half)
388 safe_deallocate_a(this%e_field_dt_full)
391 call this%get_efield(this%iteration%value(), this%e_field)
392 call this%get_efield(this%iteration%value()+algo%dt/
m_two, this%e_field_dt_half)
393 call this%get_efield(this%iteration%value()+algo%dt, this%e_field_dt_full)
398 do ip = 1, this%gr%np
400 this%e_field(ip,idim), this%gamma_p, this%omega_p, this%strength_p)
403 this%e_field_dt_half(ip, idim), this%gamma_p, this%omega_p, this%strength_p)
406 this%e_field_dt_half(ip, idim), this%gamma_p, this%omega_p, this%strength_p)
409 this%e_field_dt_full(ip, idim), this%gamma_p, this%omega_p, this%strength_p)
411 this%current_p(ip, idim) = this%current_p(ip, idim) + algo%dt / 6.0_real64 * &
412 (k1(idim) +
m_two * (k2(idim) + k3(idim)) + k4(idim))
415 updated_quantities = [
"current"]
422 call profiling_out(trim(this%namespace%get())//
":"//trim(operation%id))
427 real(real64),
intent(in) :: current_p
428 real(real64),
intent(in) :: e_field
429 real(real64),
intent(in) :: gamma_p
430 real(real64),
intent(in) :: omega_p
431 real(real64),
intent(in) :: strength_p
432 real(real64) :: current_dot
434 current_dot = - gamma_p * current_p + strength_p *
p_ep * omega_p**2 * e_field
442 real(real64),
intent(in) :: tol
459 call profiling_in(trim(partner%namespace%get())//
":"//
"COPY_QUANTITY_INTER")
461 select type (interaction)
463 interaction%partner_field(:,:) = partner%current_p
464 call interaction%do_mapping()
466 message(1) =
"Unsupported interaction."
470 call profiling_out(trim(partner%namespace%get())//
":"//
"COPY_QUANTITY_INTER")
478 character(len=256) :: filename
484 call profiling_in(trim(this%namespace%get())//
":"//
"RESTART_WRITE")
487 do idir = 1, this%space%dim
488 write(filename,
"(A,A)")
"current_p-",
index2axis(idir)
490 this%current_p(:, idir), err)
493 call iter%start(this%interactions)
494 do while (iter%has_next())
495 interaction => iter%get_next()
496 select type (interaction)
498 call interaction%write_restart(this%gr, this%space, this%restart_dump, err)
503 message(1) =
"Successfully wrote restart data for system "//trim(this%namespace%get())
509 call profiling_out(trim(this%namespace%get())//
":"//
"RESTART_WRITE")
518 character(len=256) :: filename
524 call profiling_in(trim(this%namespace%get())//
":"//
"RESTART_READ")
529 write(filename,
"(A,A)")
"current_p-",
index2axis(idir)
531 this%current_p(:, idir), err)
534 call iter%start(this%interactions)
535 do while (iter%has_next())
536 interaction => iter%get_next()
537 select type (interaction)
539 call interaction%read_restart(this%gr, this%space, this%restart_load, err)
545 this%from_scratch = .false.
548 this%current_p(:, :) =
m_zero
552 call profiling_out(trim(this%namespace%get())//
":"//
"RESTART_READ")
564 this%kinetic_energy =
m_zero
574 integer :: first, id, idir
575 character(len=130) :: aux
581 select type (algo => this%algo)
584 if (this%iteration%counter() == 0)
then
587 first = this%iteration%counter() + 1
590 call io_mkdir(
'td.general', this%namespace)
592 trim(
io_workpath(
"td.general/current_at_points", this%namespace)))
595 if (this%iteration%counter() == 0)
then
598 '################################################################################')
611 do id = 1, this%n_output_points
613 write(aux,
'(a,i1,a,i1,a)')
'j(', id,
',', idir,
')'
623 do id = 1, this%n_output_points
630 '################################################################################')
635 if (first == 0)
call this%output_write()
648 real(real64) :: dmin, dtmp(3)
649 integer :: ip, pos_index, rankmin
652 call profiling_in(trim(this%namespace%get())//
":"//
"OUTPUT_WRITE")
654 select type (algo => this%algo)
656 do ip = 1, this%n_output_points
657 pos_index =
mesh_nearest_point(this%gr, this%selected_points_coordinate(ip,:), dmin, rankmin)
658 if (this%gr%mpi_grp%rank == rankmin)
then
659 dtmp(:) = this%current_p(pos_index,:)
661 if (this%gr%parallel_in_domains)
then
662 call this%gr%mpi_grp%bcast(dtmp(:), 3, mpi_double_precision, rankmin)
668 call profiling_out(trim(this%namespace%get())//
":"//
"OUTPUT_WRITE")
674 do ip = 1, this%n_output_points
675 dtmp = this%current_at_point(ip,1:3)
684 call profiling_out(trim(this%namespace%get())//
":"//
"OUTPUT_WRITE")
707 real(real64),
intent(in) :: time
708 real(real64),
contiguous,
intent(inout) :: efield(:, :)
711 real(real64),
allocatable :: efield_tmp(:, :)
715 safe_allocate(efield_tmp(1:this%gr%np, 1:3))
718 call iter%start(this%interactions)
719 do while (iter%has_next())
720 select type (interaction => iter%get_next())
722 call interaction%interpolate(time, efield_tmp)
726 safe_deallocate_a(efield_tmp)
737 safe_deallocate_a(this%current_p)
738 safe_deallocate_a(this%e_field)
739 safe_deallocate_a(this%selected_points_coordinate)
740 safe_deallocate_a(this%current_at_point)
real(real64) function current_derivative_dir(current_p, e_field, gamma_p, omega_p, strength_p)
constant times a vector plus a vector
Writes to the corresponding file and adds one to the iteration. Must be called after write_iter_init(...
This module defines the abstract interfact for algorithm factories.
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.
subroutine dispersive_medium_init_parallelization(this, grp)
logical function dispersive_medium_is_tolerance_reached(this, tol)
logical function dispersive_medium_do_algorithmic_operation(this, operation, updated_quantities)
subroutine dispersive_medium_init_interaction(this, interaction)
subroutine dispersive_medium_restart_write_data(this)
logical function dispersive_medium_restart_read_data(this)
subroutine dispersive_medium_finalize(this)
subroutine dispersive_medium_init_interaction_as_partner(partner, interaction)
class(dispersive_medium_t) function, pointer dispersive_medium_constructor(namespace)
The factory routine (or constructor) allocates a pointer of the corresponding type and then calls the...
subroutine, public dispersive_medium_init(this, namespace)
The init routine is a module level procedure This has the advantage that different classes can have d...
subroutine dispersive_medium_output_start(this)
subroutine dispersive_medium_update_kinetic_energy(this)
subroutine dispersive_medium_initialize(this)
subroutine dispersive_medium_get_efield(this, time, efield)
subroutine dispersive_medium_output_finish(this)
subroutine dispersive_medium_output_write(this)
subroutine dispersive_medium_copy_quantities_to_interaction(partner, interaction)
This module implements the field transfer.
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public p_ep
complex(real64), parameter, public m_z0
real(real64), parameter, public m_one
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 mxll_e_field_to_matter
integer, parameter, public current_to_mxll_field
This module defines the abstract interaction_t class, and some auxiliary classes for interactions.
This module defines classes and functions for interaction partners.
character(len=max_path_len) function, public io_workpath(path, namespace)
subroutine, public io_mkdir(fname, namespace, parents)
This module defines a linear medium for use in classical electrodynamics calculations.
This module defines the meshes, which are used in Octopus.
integer function, public mesh_nearest_point(mesh, pos, dmin, rankmin)
Returns the index of the point which is nearest to a given vector position pos.
subroutine, public messages_not_implemented(feature, 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)
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
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
integer, parameter, public mxll_field_total
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.
character(len=algo_label_len), parameter, public store_current_status
character(len=algo_label_len), parameter, public rk4_finish
character(len=algo_label_len), parameter, public rk4_start
character(len=algo_label_len), parameter, public rk4_propagate
character(len=algo_label_len), parameter, public rk4_extrapolate
This module defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
Implementation details for regridding.
subroutine, public restart_init(restart, namespace, data_type, type, mc, ierr, mesh, dir, exact)
Initializes a restart object.
integer, parameter, public restart_type_dump
logical pure function, public restart_skip(restart)
Returns true if the restart information should neither be read nor written. This might happen because...
integer, parameter, public restart_td
integer, parameter, public restart_type_load
This module implements the abstract system type.
subroutine, public system_init_parallelization(this, grp)
Basic functionality: copy the MPI group. This function needs to be implemented by extended types that...
subroutine, public system_end(this)
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_out
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
This module is intended to contain simple general-purpose utility functions and procedures.
character pure function, public index2axis(idir)
Explicit interfaces to C functions, defined in write_iter_low.cc.
Descriptor of one algorithmic operation.
Class to transfer a current to a Maxwell field.
dispersive medium for classical electrodynamics calculations
These class extend the list and list iterator to make an interaction list.
abstract interaction class
surrogate interaction class to avoid circular dependencies between modules.
This is defined even when running serial.
class to transfer a Maxwell field to a medium
Abstract class implementing propagators.
Implements a the 4th order Runge Kutta propagator.
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...
Abstract class for systems.