Octopus
external_waves.F90
Go to the documentation of this file.
1!! Copyright (C) 2023 E.I. Albar, F. Bonafe and Heiko Appel
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17
18#include "global.h"
19
21 use accel_oct_m
24 use clock_oct_m
26 use debug_oct_m
30 use global_oct_m
31 use grid_oct_m
32 use index_oct_m
37 use, intrinsic :: iso_fortran_env
39 use io_oct_m
41 use lasers_oct_m
45 use mesh_oct_m
47 use mpi_oct_m
52 use parser_oct_m
56 use string_oct_m
57 use types_oct_m
58 use unit_oct_m
62
63 implicit none
64
65 private
66 public :: &
73
74 type bessel_beam_t
75 integer, allocatable :: helicity(:)
76 integer, allocatable :: m_order(:)
77 real(real64), allocatable :: amp(:)
78 real(real64), allocatable :: theta_k(:)
79 real(real64), allocatable :: omega(:)
80 real(real64), allocatable :: shift(:,:)
81 logical, allocatable :: envelope(:)
82 integer, allocatable :: lin_dir(:)
83 contains
84 procedure :: init => bessel_beam_init
85 procedure :: function => bessel_beam_function
87 end type bessel_beam_t
88
90 integer :: points_number
91 integer, allocatable :: points_map(:)
92 integer :: number
93 integer, allocatable :: modus(:)
94 integer, allocatable :: field_type(:)
95 character(len=1024), allocatable :: e_field_string(:,:)
96 real(real64), allocatable :: k_vector(:,:)
97 real(real64), allocatable :: v_vector(:,:)
98 complex(real64), allocatable :: e_field(:,:)
99 real(real64), allocatable :: pw_phase(:)
100 type(mxf_t), allocatable :: mx_function(:)
101 integer :: out_file
102 logical :: output_from_point = .false.
103 real(real64), allocatable :: selected_point_coordinate(:)
104 real(real64), allocatable :: selected_point_field(:)
105 real(real64) :: c_factor
106 type(accel_mem_t) :: buff_map
107 type(bessel_beam_t) :: bessel
108 contains
109 procedure :: init_interaction_as_partner => external_waves_init_interaction_as_partner
110 procedure :: update_quantity => external_waves_update_quantity
111 procedure :: copy_quantities_to_interaction => external_waves_copy_quantities_to_interaction
112 final :: external_waves_end
114
115 interface external_waves_t
116 module procedure external_waves_constructor
117 end interface external_waves_t
118
119contains
120
121 function external_waves_constructor(namespace) result(this)
122 class(external_waves_t), pointer :: this
123 type(namespace_t), intent(in) :: namespace
124
125 integer :: iq
126 character(len=:), allocatable :: quantities(:)
127
129
130 safe_allocate(this)
131
132 this%namespace = namespace_t("ExternalSource", parent=namespace)
133
134 message(1) = 'Plane-wave is currently always 3D and non-periodic.'
135 call messages_warning(1)
136 call external_waves_init(this, this%namespace)
137
138 quantities = [character(16) :: "E field", "vector potential", "B field"]
139 do iq = 1, size(quantities)
140 call this%quantities%add(quantity_t(quantities(iq), always_available = .true., &
141 updated_on_demand = .true., iteration = clock_t()))
142 end do
143
144 this%supported_interactions_as_partner = [mxll_e_field_to_matter, mxll_b_field_to_matter, &
146
148 end function external_waves_constructor
149
150 ! ---------------------------------------------------------
151 subroutine external_waves_init_interaction_as_partner(partner, interaction)
152 class(external_waves_t), intent(in) :: partner
153 class(interaction_surrogate_t), intent(inout) :: interaction
154
156
157 select type (interaction)
158 type is (lorentz_force_t)
159 ! Nothing to be initialized
161 ! Nothing to be initialized
163 ! Nothing to be initialized
165 ! Nothing to be initialized
166 class default
167 message(1) = "Unsupported interaction."
168 call messages_fatal(1, namespace=partner%namespace)
169 end select
173 ! ---------------------------------------------------------
174 subroutine external_waves_update_quantity(this, label)
175 class(external_waves_t), intent(inout) :: this
176 character(len=*), intent(in) :: label
180 select case (label)
181 case ("E field", "B field", "vector potential")
182 ! We will not update the quantities here because they are computed
183 ! on-the-fly when copying them to the corresponding interaction (see
184 ! copy_quantities_to_interaction routine)
185 case default
186 message(1) = "Incompatible quantity."
187 call messages_fatal(1, namespace=this%namespace)
188 end select
193 ! ---------------------------------------------------------
194 subroutine external_waves_copy_quantities_to_interaction(partner, interaction)
195 class(external_waves_t), intent(inout) :: partner
196 class(interaction_surrogate_t), intent(inout) :: interaction
198 class(quantity_t), pointer :: quantity
201
202 select type(interaction)
204 quantity => partner%quantities%get("E field")
205 interaction%system_field = m_zero
206 call external_waves_eval(partner, quantity%iteration%value(), interaction%system_gr, "E field", &
207 interaction%system_field)
208 call interaction%do_mapping()
209
210 class is (mxll_vec_pot_to_matter_t)
211 quantity => partner%quantities%get("vector potential")
212 interaction%system_field = m_zero
213 call external_waves_eval(partner, quantity%iteration%value(), interaction%system_gr, &
214 "vector potential", interaction%system_field)
215 call interaction%do_mapping()
216
217 class is (mxll_b_field_to_matter_t)
218 quantity => partner%quantities%get("B field")
219 interaction%system_field = m_zero
220 call external_waves_eval(partner, quantity%iteration%value(), interaction%system_gr, "B field", &
221 interaction%system_field, der=interaction%system_gr%der)
222 call interaction%do_mapping()
223
224 class default
225 message(1) = "Incompatible interaction."
226 call messages_fatal(1, namespace=partner%namespace)
227 end select
228
231
232 ! ---------------------------------------------------------
233 ! Load the external source for the multisystem framework
234 subroutine load_external_waves(partners, namespace)
235 class(partner_list_t), intent(inout) :: partners
236 type(namespace_t), intent(in) :: namespace
237
238 logical :: has_source
239
240 push_sub(load_external_waves)
241
242 !%Variable AnalyticalExternalSource
243 !%Type logical
244 !%Default no
245 !%Section Maxwell
246 !%Description
247 !% This means the analytical evaluation of formula will be used, Maxwell propagation will not be used.
248 !%End
249 call parse_variable(namespace, 'AnalyticalExternalSource', .false., has_source)
250
251 if (has_source) then
252 call partners%add(external_waves_t(namespace))
253 end if
254
255 pop_sub(load_external_waves)
256 end subroutine load_external_waves
257
258
260 ! ---------------------------------------------------------
261 subroutine external_waves_init(external_waves, namespace)
262 type(external_waves_t), intent(inout) :: external_waves
263 type(namespace_t), intent(in) :: namespace
264 type(block_t) :: blk
265 integer :: il, nlines, ncols, iex_norm, idim
266 integer, parameter :: sys_dim = 3
267 real(real64) :: k_vector(sys_dim), velocity(sys_dim), x_pos(sys_dim)
268 real(real64) :: x_norm, dummy(sys_dim), k_dot_e , test_limit, k_norm, output_pos(3)
269 complex(real64) :: e_field(sys_dim)
270 character(len=1024) :: k_string(sys_dim)
271 character(len=1), dimension(sys_dim), parameter :: dims = ["x", "y", "z"]
272 character(len=1024) :: mxf_expression
273
274 push_sub(external_waves_init)
275
276 call profiling_in('EXTERNAL_WAVES_INIT')
277
278 test_limit = 10.0e-9_real64
279
280 !%Variable ExternalSourceBesselOutput
281 !%Type block
282 !%Section Maxwell
283 !%Description
284 !% The ExternalSourceBesselOutput block allows to output analytically calculated fields at a
285 !% particular point in space. The columns denote the x, y, and z coordinate of the point.
286 !% Please be aware that ExternalSource lives on the grid of the system that it is applied to.
287 !% Therefore, it might not be evaluated at every point in space. When comparing, please be sure
288 !% to check the log and compare if your required point in space matches the evaluated position.
289 !%
290 !% <tt>%ExternalSourceBesselOutput
291 !% <br>&nbsp;&nbsp; -1.0 | 2.0 | 4.0
292 !% <br>%</tt>
293 !%
294 !%End
295
296 if (parse_block(namespace, 'ExternalSourceBesselOutput', blk) == 0) then
297 nlines = parse_block_n(blk)
298 if (nlines > 1 ) then
299 message(2) = 'ExternalSource output is limited to one point.'
300 call messages_fatal(1, namespace=namespace)
301 end if
302 ncols = parse_block_cols(blk,0)
303 if (ncols /= 3 ) then
304 message(1) = 'ExternalSourceBesselOutput must have 3 columns.'
305 call messages_fatal(1, namespace=namespace)
306 end if
307 external_waves%output_from_point= .true.
308 safe_allocate(external_waves%selected_point_coordinate(1:3))
309 safe_allocate(external_waves%selected_point_field(1:3))
310
311 do idim = 1, 3
312 call parse_block_float(blk, 0, idim-1, output_pos(idim), units_inp%length)
313 end do
314 external_waves%selected_point_coordinate(1:3) = output_pos(1:3)
315 external_waves%selected_point_field(1:3) = m_zero
316
317 call parse_block_end(blk)
318 call io_mkdir('ExternalSource')
319 external_waves%out_file = io_open('bessel_source_at_point', namespace=namespace, action='write')
320 write(external_waves%out_file, '(12a) ') '# Time (a.u.) ' , ' Field-x ' , ' Field-y ' , ' Field-z '
321
322 else
323 external_waves%output_from_point= .false.
324
325 end if
326
327 ! This variable is documented in hamiltonian_mxll.F90
328 call parse_variable(namespace, 'SpeedOfLightFactor', m_one, external_waves%c_factor)
329
330 !%Variable MaxwellIncidentWaves
331 !%Type block
332 !%Section Maxwell
333 !%Description
334 !% The initial electromagnetic fields can be set by the user
335 !% with the <tt>MaxwellIncidentWaves</tt> block variable.
336 !% The electromagnetic fields have to fulfill the
337 !% Maxwells equations in vacuum.
338 !% For a Maxwell propagation, setting the electric field is sufficient,
339 !% the magnetic field (for plane waves) will be calculated from it as 1/(c.|k|) . (k x E).
340 !%
341 !% Example:
342 !%
343 !% <tt>%MaxwellIncidentWaves
344 !% <br>&nbsp;&nbsp; plane_wave_parser | "field_type"| "k1x" | "k1y" | "k1z" | "E1x" | "E1z" | "E1x"
345 !% <br>&nbsp;&nbsp; plane_wave_mx_function | "field_type"| "E4x" | "E4y" | "E4z" | mx_envelope_name | phase
346 !% <br>&nbsp;&nbsp; bessel_function | "field_type"| A_0 | m | omega | helicity | <math>\theta_{k}</math> | mx_envelope_name | lin_dir
347 !% <br>%</tt>
348 !%
349 !% Field type can be "electric_field" or "vector_potential". Note that in order to couple to an
350 !% electronic system, the <tt>MaxwellCouplingMode</tt> variable needs to be set to a coupling type
351 !% compatible with the requested field type ("electric_field" is compatible with length gauge,
352 !% while "vector_potential" is compatible with velocity gauge and full minimal coupling).
353 !% Otherwise, the field will not be calculated or applied to the electronic Hamiltonian.
354 !%
355 !%Option plane_wave_parser 0
356 !% Parser input modus
357 !%Option plane_wave_mx_function 1
358 !% The incident wave envelope is defined by an mx_function
359 !%Option bessel_function 2
360 !% The incident source is a generalized Bessel beam, parametrized by its amplitude, opening angle, helicity, order and frequency.
361 !% This beam is a solution of Maxwell equations, and inherently circularly polarized and is parametrized by its amplitude,
362 !% opening angle, helicity, order and frequency.
363 !% Please keep in mind, if you set linear polarization lin_dir,
364 !% you will obtain a linearly polarized Bessel beam.
365 !%End
366
367 if (parse_block(namespace, 'MaxwellIncidentWaves', blk) == 0) then
368
369 call messages_print_with_emphasis(msg='Substitution of the electromagnetic incident waves', namespace=namespace)
370
371 ! find out how many lines (i.e. states) the block has
372 nlines = parse_block_n(blk)
373
374 external_waves%number = nlines
375 safe_allocate(external_waves%modus(1:nlines))
376 safe_allocate(external_waves%e_field_string(1:3, 1:nlines))
377 safe_allocate(external_waves%e_field(1:3, 1:nlines))
378 safe_allocate(external_waves%k_vector(1:3, 1:nlines))
379 safe_allocate(external_waves%v_vector(1:3, 1:nlines))
380 safe_allocate(external_waves%mx_function(1:nlines))
381 safe_allocate(external_waves%field_type(1:nlines))
382 safe_allocate(external_waves%pw_phase(1:nlines))
383 external_waves%pw_phase = m_zero
384
385 call external_waves%bessel%init(nlines, 3)
386
387 do il = 1, nlines
388 ncols = parse_block_cols(blk, il - 1)
389 if ((ncols < 5) .or. (ncols > 9)) then
390 message(1) = 'Each line in the MaxwellIncidentWaves block must have five to nine columns.'
391 call messages_fatal(1, namespace=namespace)
392 end if
393
394 ! check input modus e.g. parser of defined functions
395 call parse_block_integer(blk, il - 1, 0, external_waves%modus(il))
396 call parse_block_integer(blk, il - 1, 1, external_waves%field_type(il))
397
398 ! parse formula string
399 if (external_waves%modus(il) == option__maxwellincidentwaves__plane_wave_parser) then
400 do idim = 1, 3
401 call parse_block_string( blk, il - 1, idim + 1, k_string(idim))
402 call parse_block_string( blk, il - 1, 3 + idim + 1, external_waves%e_field_string(idim, il))
403 end do
404 write(message(1), '(a,i2,a) ') 'Substituting electromagnetic incident wave ', il, ' with the expressions: '
405 call messages_info(1, namespace=namespace)
406 do idim = 1, 3
407 write(message(idim), '(6a)') ' Wave vector k('//dims(idim)//') = ', trim(k_string(idim))
408 write(message(idim+1), '(2a)') ' E-field('//dims(idim)//') for t_0 = ', trim(external_waves%e_field_string(idim, il))
409 end do
410 call messages_info(6, namespace=namespace)
411
412 do idim = 1, 3
413 call conv_to_c_string(k_string(idim))
414 call conv_to_c_string(external_waves%e_field_string(idim, il))
415 end do
416
417 x_pos(:) = m_zero
418 x_norm = m_zero
419 do idim = 1, 3
420 call parse_expression(k_vector(idim), dummy(idim), idim, x_pos, x_norm, m_zero, k_string(idim))
421 end do
422
423 k_vector = units_to_atomic(unit_one/units_inp%length, k_vector)
424 k_norm = norm2(k_vector)
425
426 velocity(:) = k_vector(:) / k_norm * p_c * external_waves%c_factor
427 external_waves%k_vector(:,il) = k_vector(:)
428 external_waves%v_vector(:,il) = velocity(:)
429
430 else if (external_waves%modus(il) == option__maxwellincidentwaves__plane_wave_mx_function) then
431 do idim = 1, 3
432 call parse_block_cmplx( blk, il - 1, idim + 1, e_field(idim))
433 end do
434 call parse_block_string( blk, il - 1, 3 + 2, mxf_expression)
435
436 write(message(1), '(a,i2) ') 'Substituting electromagnetic incident wave ', il
437 write(message(2), '(a)' ) 'with the expression: '
438 call messages_info(2, namespace=namespace)
439
440 do idim = 1, 3
441 write(message(idim), '(a,f9.4,sp,f9.4,"i")') ' E-field('//trim(dims(idim))//') complex amplitude = ', &
442 real(e_field(idim)), aimag(e_field(idim))
443 end do
444 write(message(4), '(2a)' ) ' Maxwell wave function name = ', trim(mxf_expression)
445 call messages_info(4, namespace=namespace)
446
447 call mxf_read(external_waves%mx_function(il), namespace, trim(mxf_expression), iex_norm)
448 if (iex_norm /= 0) then
449 write(message(1),'(3A)') 'Ex_norm in the ""', trim(mxf_expression), &
450 '"" field defined in the MaxwellIncidentWaves block'
451 call messages_fatal(1, namespace=namespace)
452 end if
453 if (parse_block_cols(blk, il-1) == 7) then
454 call parse_block_float( blk, il - 1, 3 + 3 , external_waves%pw_phase(il))
455 end if
456 e_field = units_to_atomic(units_inp%energy/units_inp%length, e_field)
457 k_vector(1:3) = external_waves%mx_function(il)%k_vector(1:3)
458 k_norm = norm2(k_vector)
459
460 k_dot_e = real(dot_product(k_vector, e_field), real64)
461 if (abs(k_dot_e) > test_limit) then
462 write(message(1), '(a) ') 'The wave vector k or its electric field are not perpendicular. '
463 write(message(2), '(a,f8.3,a)' ) 'Their dot product yields the magnitude', abs(k_dot_e) , ' while'
464 write(message(3), '(a,f8.3,a)' ) 'tolerance is ', test_limit ,'.'
465 call messages_fatal(3, namespace=namespace)
466 end if
467 if (k_norm < 1e-10) then
468 message(1) = 'The k vector is not set correctly: |k|~0 .'
469 call messages_fatal(1, namespace=namespace)
470 end if
471
472 external_waves%e_field(:,il) = e_field(:)
473 external_waves%k_vector(:,il) = k_vector(:)
474 external_waves%v_vector(:,il) = k_vector(:) / k_norm * p_c * external_waves%c_factor
475
476 else if (external_waves%modus(il) == option__maxwellincidentwaves__bessel_function) then
477 call parse_block_float( blk, il - 1, 2 , external_waves%bessel%amp(il))
478 call parse_block_integer( blk, il - 1, 3 , external_waves%bessel%m_order(il))
479 call parse_block_float( blk, il - 1, 4 , external_waves%bessel%omega(il))
480 call parse_block_integer( blk, il - 1, 5 , external_waves%bessel%helicity(il))
481 call parse_block_float( blk, il - 1, 6 , external_waves%bessel%theta_k(il))
482 if (parse_block_cols(blk, il-1) > 7) then
483 call parse_block_string( blk, il - 1, 7 , mxf_expression)
484 external_waves%bessel%envelope(il) = .true.
485 call mxf_read(external_waves%mx_function(il), namespace, trim(mxf_expression), iex_norm)
486 end if
487 if (parse_block_cols(blk, il-1) == 9) then
488 call parse_block_integer( blk, il - 1, 8 , external_waves%bessel%lin_dir(il))
489 end if
490
491 write(message(1), '(a,i2) ') 'Incident Bessel Beam', il
492 call messages_info(1, namespace=namespace)
493
494 if (abs(external_waves%bessel%helicity(il)) /= 1) then
495 write(message(1),'(A)') 'Helicity has to be either +1 or -1 !'
496 call messages_fatal(1, namespace=namespace)
497 end if
498
499 write(message(1), '(a,f5.3)' ) ' Bessel Amplitude ', external_waves%bessel%amp(il)
500 write(message(2), '(a,i2)' ) ' Bessel Order m', external_waves%bessel%m_order(il)
501 write(message(3), '(a,f5.3)' ) ' Bessel Frequency ', external_waves%bessel%omega(il)
502 write(message(4), '(a,i2)' ) ' Bessel Helicity ', external_waves%bessel%helicity(il)
503 write(message(5), '(a,f5.3)' ) ' Bessel Opening Angle ', external_waves%bessel%theta_k(il)
504 call messages_info(4, namespace=namespace)
505
506 if (external_waves%bessel%lin_dir(il)/= 0) then
507 write(message(5), '(a,i2)' ) ' Bessel is Linearly Polarized in Direction : ', external_waves%bessel%lin_dir(il)
508 call messages_info(4, namespace=namespace)
509 end if
510
511 end if
512 end do
513
514 call parse_block_end(blk)
515
516 call messages_print_with_emphasis(namespace=namespace)
517 else
518 external_waves%number = 0
519
520 end if
521
522 !%Variable BesselBeamAxisShift
523 !%Type block
524 !%Section Maxwell
525 !%Description
526 !% The BesselBeamAxisShift block allows to shift the Bessel Beam, which is centered at (0,0,0) as default.
527 !% Selected position point will be used as the new center of the Bessel Beam.
528 !% When defining a BesselBeamAxisShift, please make sure to define a shift for each Bessel source you use,
529 !% then it is possible to tell which source is shifted according to which BesselShift, respectively.
530 !% <tt>%BesselBeamAxisShift
531 !% <br>&nbsp;&nbsp; 0.0 | 2.0 | 5.0
532 !% <br>%</tt>
533 !%
534 !%End
535
536 if (parse_block(namespace, 'BesselBeamAxisShift', blk) == 0) then
537 nlines = parse_block_n(blk)
538 ncols = parse_block_cols(blk,0)
539 if (ncols /= 3 ) then
540 message(1) = 'BesselBeamAxisShift must have 3 columns.'
541 call messages_fatal(1, namespace=namespace)
542 end if
543
544 do il = 1, nlines
545 do idim = 1, 3
546 call parse_block_float(blk, 0, idim-1, external_waves%bessel%shift(il, idim), units_inp%length)
547 end do
548 end do
549
550 call parse_block_end(blk)
551 end if
552
553 call profiling_out('EXTERNAL_WAVES_INIT')
554
555 pop_sub(external_waves_init)
556 end subroutine external_waves_init
557
558 ! ---------------------------------------------------------
559 subroutine external_waves_end(external_waves)
560 type(external_waves_t), intent(inout) :: external_waves
561
562 push_sub(external_waves_end)
563
564 if (external_waves%output_from_point) then
565 call io_close(external_waves%out_file)
566 safe_deallocate_a(external_waves%selected_point_coordinate)
567 safe_deallocate_a(external_waves%selected_point_field)
568 end if
569
570 safe_deallocate_a(external_waves%bessel%shift)
571 safe_deallocate_a(external_waves%points_map)
572 safe_deallocate_a(external_waves%modus)
573 safe_deallocate_a(external_waves%e_field_string)
574 safe_deallocate_a(external_waves%k_vector)
575 safe_deallocate_a(external_waves%v_vector)
576 safe_deallocate_a(external_waves%e_field)
577 safe_deallocate_a(external_waves%mx_function)
578 safe_deallocate_a(external_waves%pw_phase)
579
580 if (accel_is_enabled()) then
581 call accel_release_buffer(external_waves%buff_map)
582 end if
583
584 pop_sub(external_waves_end)
585 end subroutine external_waves_end
586
587 ! ---------------------------------------------------------
589 subroutine external_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
590 type(external_waves_t), intent(inout) :: external_waves
591 real(real64), intent(in) :: time
592 class(mesh_t), intent(in) :: mesh
593 character(len=*), intent(in) :: type_of_field
594 real(real64), intent(out) :: out_field_total(:, :)
595 type(derivatives_t), optional, intent(in):: der
596
597
598 push_sub(external_waves_eval)
599
600 call profiling_in('EXTERNAL_WAVES_EVAL')
601
602 out_field_total = m_zero
603
604 call plane_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der=der)
605 call bessel_source_eval(external_waves, time, mesh, type_of_field, out_field_total, der=der)
606
607 call profiling_out('EXTERNAL_WAVES_EVAL')
608
609 pop_sub(external_waves_eval)
610 end subroutine external_waves_eval
611
612 ! ---------------------------------------------------------
614 subroutine plane_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
615 type(external_waves_t), intent(inout) :: external_waves
616 real(real64), intent(in) :: time
617 class(mesh_t), intent(in) :: mesh
618 character(len=*), intent(in) :: type_of_field
619 real(real64), intent(out) :: out_field_total(:, :)
620 type(derivatives_t), optional, intent(in):: der
621
622 integer :: wn
623 real(real64), allocatable :: pw_field(:,:), ztmp(:,:), b_field_aux(:,:)
624 real(real64) :: p_c_
625 integer, allocatable :: indices_pw_parser(:)
626 integer, allocatable :: indices_mx_ftc(:)
627 integer :: n_plane_waves, n_points
628
629 push_sub(plane_waves_eval)
630
631 call profiling_in('PLANE_WAVES_EVAL')
632
633 indices_pw_parser = pack([(wn, wn = 1,external_waves%number)], &
634 external_waves%modus == option__maxwellincidentwaves__plane_wave_parser)
635
636 indices_mx_ftc = pack([(wn, wn = 1,external_waves%number)], &
637 external_waves%modus == option__maxwellincidentwaves__plane_wave_mx_function)
638
639 n_plane_waves = size(indices_pw_parser) + size(indices_mx_ftc)
640
641 p_c_ = p_c * external_waves%c_factor
642
643 if (n_plane_waves == 0) then
644 call profiling_out('PLANE_WAVES_EVAL')
645 pop_sub(plane_waves_eval)
646 return
647 end if
648
649 if (type_of_field == "B field" .and. any(external_waves%field_type == e_field_vector_potential)) then
650 assert(present(der))
651 safe_allocate(ztmp(mesh%np, size(out_field_total, dim=2)))
652 n_points = mesh%np_part
653 else
654 n_points = mesh%np
655 end if
656 safe_allocate(pw_field(n_points, size(out_field_total, dim=2)))
657 pw_field(:,:) = m_zero
658
659 ! The E_field (or A_field, rescaled later) we calculate always
660 do wn = 1, external_waves%number
661
662 select case(external_waves%modus(wn))
663 case (option__maxwellincidentwaves__plane_wave_parser)
664 call pw_parsed_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
665
666 case (option__maxwellincidentwaves__plane_wave_mx_function)
667 call pw_mx_function_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
668 end select
669
670 select case (external_waves%field_type(wn))
671
672 case(e_field_electric)
673
674 select case (type_of_field)
675 case ("E field")
676 out_field_total(1:mesh%np,:) = out_field_total(1:mesh%np,:) + pw_field(1:mesh%np,:)
677 case ("vector potential")
678 call messages_not_implemented("Calculation of a vector potential from a plane wave specified as electric field")
679 case ("B field")
680 safe_allocate(b_field_aux(1:mesh%np, 1:mesh%box%dim))
681 call get_pw_b_field(external_waves, mesh, wn, pw_field, b_field_aux)
682 out_field_total(:,:) = out_field_total(:,:) + b_field_aux(:,:)
683 safe_deallocate_a(b_field_aux)
684 end select
685
687
688 select case (type_of_field)
689 case ("E field")
690 call messages_not_implemented("Calculation of an electric field from a plane wave specified as vector potential")
691 case ("vector potential")
692 out_field_total(1:mesh%np,:) = out_field_total(1:mesh%np,:) - m_one/p_c_ * pw_field(1:mesh%np,1:3)
693 case ("B field")
694 call dderivatives_curl(der, pw_field(1:mesh%np_part,1:3), ztmp(1:mesh%np,1:3), set_bc = .false.)
695 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c_ * ztmp(1:mesh%np, 1:3)
696 end select
697
698 end select
699 end do
700
701 safe_deallocate_a(pw_field)
702 safe_deallocate_a(ztmp)
703 call profiling_out('PLANE_WAVES_EVAL')
704
705 pop_sub(plane_waves_eval)
706
707 end subroutine plane_waves_eval
708
709 ! ---------------------------------------------------------
711 subroutine pw_parsed_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
712 type(external_waves_t), intent(inout) :: external_waves
713 integer, intent(in) :: wn
714 real(real64), intent(in) :: time
715 class(mesh_t), intent(in) :: mesh
716 integer, intent(in) :: n_points
717 real(real64), intent(out) :: pw_field(:,:)
718
719 real(real64) :: x_prop(3), x_norm
720 real(real64) :: velocity_time(3)
721 real(real64) :: parsed_field(3)
722 real(real64) :: dummy(3)
723 integer :: idim, ip
724
725 velocity_time(:) = external_waves%v_vector(1:3, wn) * time
726 do idim = 1, 3
727 call parse_expression(parsed_field(idim), dummy(idim), 3, x_prop, x_norm, m_zero, &
728 external_waves%e_field_string(idim, wn))
729 do ip = 1, n_points
730 x_prop = mesh%x(ip, :) - velocity_time
731 x_norm = norm2(x_prop(1:3))
732 pw_field(ip, idim) = units_to_atomic(units_inp%energy/units_inp%length, parsed_field(idim))
733 end do
734 end do
735
736 end subroutine pw_parsed_evaluation
737
738 ! ---------------------------------------------------------
740 subroutine pw_mx_function_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
741 type(external_waves_t), intent(inout) :: external_waves
742 integer, intent(in) :: wn
743 real(real64), intent(in) :: time
744 class(mesh_t), intent(in) :: mesh
745 integer, intent(in) :: n_points
746 real(real64), intent(out) :: pw_field(:,:)
747
748 real(real64) :: x_prop(3), x_norm
749 real(real64) :: velocity_time(3)
750 complex(real64) :: efield_ip(3)
751 complex(real64) :: e0(3)
752 integer :: ip
753
754 velocity_time(:) = external_waves%v_vector(1:3, wn) * time
755 e0(:) = external_waves%e_field(1:3, wn)
756 do ip = 1, n_points
757 x_prop = mesh%x(ip, :) - velocity_time
758 x_norm = norm2(x_prop(1:3))
759 efield_ip = mxf(external_waves%mx_function(wn), x_prop, external_waves%pw_phase(wn))
760 pw_field(ip, :) = real(e0(1:3) * efield_ip, real64)
761 end do
762
763 end subroutine pw_mx_function_evaluation
764
765 ! ---------------------------------------------------------
767 subroutine get_pw_b_field(external_waves, mesh, pwidx, e_field, b_field)
768 type(external_waves_t), intent(in) :: external_waves
769 class(mesh_t), intent(in) :: mesh
770 real(real64), intent(in) :: e_field(:,:)
771 real(real64), intent(out) :: b_field(:,:)
772 integer, intent(in) :: pwidx
773
774 real(real64) :: k_vector(3), k_vector_abs
775 real(real64) :: velocity(3)
776 real(real64) :: P_c_
777 complex(real64) :: e0(3)
778 integer :: ip
779
780 velocity = external_waves%v_vector(1:3, pwidx)
781 k_vector = external_waves%k_vector(1:3, pwidx)
782 k_vector_abs = norm2(k_vector(1:3))
783 e0 = external_waves%e_field(1:3, pwidx)
784 p_c_ = p_c * external_waves%c_factor
785
786 b_field = m_zero
787 do ip = 1, mesh%np
788 b_field(ip, :) = m_one/(p_c_ * k_vector_abs) * dcross_product(k_vector, e_field(ip, :))
789 end do
790
791 end subroutine get_pw_b_field
792
793 ! ---------------------------------------------------------
795 subroutine bessel_source_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
796 type(external_waves_t), intent(inout) :: external_waves
797 real(real64), intent(in) :: time
798 class(mesh_t), intent(in) :: mesh
799 character(len=*), intent(in) :: type_of_field
800 real(real64), intent(out) :: out_field_total(:, :)
801 type(derivatives_t), optional, intent(in):: der
802
803 real(real64) :: dmin, omega, k_vector(3), c_factor
804 integer :: iline, wn, pos_index, n_points, rankmin
805 real(real64), allocatable :: shift(:,:)
806 complex(real64), allocatable :: bessel_field_total(:,:), ztmp(:,:), vec_pot(:,:)
807 integer, allocatable :: indices_bessel_ftc(:)
808 type(mxf_t) :: envelope_mxf
809
810 push_sub(bessel_source_eval)
811
812 call profiling_in('BESSEL_SOURCE_EVAL')
813
814 indices_bessel_ftc = pack([(wn, wn = 1,external_waves%number)], &
815 external_waves%modus == option__maxwellincidentwaves__bessel_function)
816
817 if (size(indices_bessel_ftc) == 0) then
818 call profiling_out('BESSEL_SOURCE_EVAL')
819 pop_sub(bessel_source_eval)
820 return
821 end if
822
823 ! Check if the BesselBeamAxisShift is defined for every incoming Bessel Beam.
824 if (allocated(external_waves%bessel%shift) .and. &
825 size(external_waves%bessel%shift(:,1)) /= size(indices_bessel_ftc)) then
826 message(1) = 'Number of BesselBeamAxisShift defined in input file'
827 message(2) = 'does not match the number of Bessel beams.'
828 call messages_fatal(2)
829 end if
830
831 safe_allocate(shift(size(indices_bessel_ftc), 3))
832 if (allocated(external_waves%bessel%shift)) then
833 shift = external_waves%bessel%shift
834 else
835 shift = m_zero
836 end if
837
838 if (type_of_field == "B field") then
839 assert(present(der))
840 safe_allocate(vec_pot(mesh%np_part, size(out_field_total, dim=2)))
841 safe_allocate(ztmp(size(out_field_total, dim=1), size(out_field_total, dim=2)))
842 n_points = mesh%np_part ! needed for curl
843 else
844 n_points = mesh%np
845 end if
846
847
848 safe_allocate(bessel_field_total(1:n_points, 1:3))
849 bessel_field_total = m_zero
850
851 do iline = 1, size(indices_bessel_ftc)
852 wn = indices_bessel_ftc(iline)
853 omega = external_waves%bessel%omega(wn)
854 k_vector = external_waves%mx_function(wn)%k_vector
855 c_factor = external_waves%c_factor
856 envelope_mxf = external_waves%mx_function(wn)
857
858 call external_waves%bessel%function(wn, shift, mesh, n_points, time, k_vector, c_factor, envelope_mxf, bessel_field_total)
859
860 select case (external_waves%field_type(wn))
861
863 ! interpreting bessel_field_total as a vector potential (as requested by the user)
864 select case (type_of_field)
865 case ("E field")
866 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) + real(m_zi*omega*bessel_field_total(1:mesh%np,1:3))
867 case ("vector potential")
868 ! For the vector potential, we multiply by -1/c becuase of the electronic Hamiltonian
869 ! being in Gaussian units
870 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c * real(bessel_field_total(1:mesh%np,1:3))
871 case ("B field")
872 call zderivatives_curl(der, bessel_field_total(1:mesh%np_part,1:3), ztmp(1:mesh%np,1:3), set_bc = .false.)
873 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c * real(ztmp(1:mesh%np, 1:3))
874 end select
875
876 case(e_field_electric)
877 ! interpreting bessel_field_total as an electric field (as requested by the user)
878 select case (type_of_field)
879 case ("E field")
880 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) + real(bessel_field_total(1:mesh%np,1:3))
881 case ("vector potential")
882 ! We calculate the vector potential as real(E/i*omega),
883 ! and convert it to the proper units by multiplying by -1/c
884 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - m_one/p_c * &
885 real(bessel_field_total(1:mesh%np,1:3)/M_zI/omega)
886 case ("B field")
887 vec_pot(1:mesh%np_part,1:3) = - m_one/p_c * real(bessel_field_total(1:mesh%np_part,1:3)/m_zi/omega)
888 call zderivatives_curl(der, vec_pot(1:mesh%np_part,1:3), ztmp(1:mesh%np,1:3), set_bc = .false.)
889 out_field_total(1:mesh%np,1:3) = out_field_total(1:mesh%np,1:3) - real(ztmp(1:mesh%np, 1:3))
890 end select
891
892 end select
893 end do
894
895 if (external_waves%output_from_point) then
896 pos_index = mesh_nearest_point(mesh, external_waves%selected_point_coordinate(1:3), dmin, rankmin)
897 if (mesh%mpi_grp%rank == rankmin) then
898 external_waves%selected_point_field(:) = out_field_total(pos_index,:)
899 write(external_waves%out_file, "(4F14.8, 4x)") time, external_waves%selected_point_field(:)
900 end if
901 end if
902
903 safe_deallocate_a(shift)
904 safe_deallocate_a(ztmp)
905 safe_deallocate_a(vec_pot)
906 safe_deallocate_a(bessel_field_total)
907 call profiling_out('BESSEL_SOURCE_EVAL')
908
909 pop_sub(bessel_source_eval)
910
911 end subroutine bessel_source_eval
912
913 ! ---------------------------------------------------------
915 subroutine bessel_beam_function(this, iline, shift, mesh, n_points, time, k_vector, c_factor, envelope_mxf, bessel_field)
916 class(bessel_beam_t) :: this
917 integer, intent(in) :: iline
918 real(real64), intent(in) :: shift(:,:), time, k_vector(3), c_factor
919 class(mesh_t), intent(in) :: mesh
920 integer, intent(in) :: n_points
921 type(mxf_t), intent(in) :: envelope_mxf
922 complex(real64), intent(out) :: bessel_field(:,:)
923
924 real(real64) :: pos(3), temp, temp2, temp3, rho, phi_rho, wigner(3)
925 real(real64) :: hel, theta, omega, amp, kappa, proj, k_norm, velocity_time(3), x_prop(3)
926 complex(real64) :: efield_ip(3)
927 real(real64) :: bessel_plus, bessel_minus
928 integer :: ip, mm, pol
929
930 assert(iline <= size(this%omega))
931 hel = real(this%helicity(iline), real64)
932 theta = this%theta_k(iline)
933 mm = this%m_order(iline)
934 amp = this%amp(iline) / sqrt(m_two)
935 omega = this%omega(iline)
936 proj = omega * cos(theta) / p_c ! k_z
937 kappa = sqrt(omega**2 - (proj* p_c)**2) ! parse omega
938 ! Set Wigner Coefficients from theta
939 wigner(1) = hel * sin(theta) / sqrt(m_two) ! mu = 0
940 wigner(2) = 0.5 * (1 + hel * cos(theta)) ! mu = 1
941 wigner(3) = 0.5 * (1 - hel * cos(theta)) ! mu = -1
942 proj = omega * cos(theta) / p_c ! k_z
943 pol = this%lin_dir(iline) ! Incoming polarization corresponding to beam in question
944
945 do ip = 1, n_points
946 pos(:) = mesh%x(ip, :) - shift(iline,:)
947 rho = norm2(pos(1:2))
948 phi_rho = atan2(pos(2) , pos(1))
949 temp = proj * pos(3) + phi_rho * (mm + 1) - omega*time ! temp, temp2 and temp3 should be unitless
950 temp2 = proj * pos(3) + phi_rho * (mm - 1) - omega*time
951 temp3 = proj * pos(3) + phi_rho * mm - omega*time
952 bessel_plus = loct_bessel(mm+1, kappa * rho / p_c)
953 bessel_minus = loct_bessel(mm-1, kappa * rho / p_c)
954
955 ! Calculate complex Ax component, if generalized bessel OR x -polarized bessel
956 if (pol /= 2) then
957 bessel_field(ip, 1) = amp * (exp(m_zi*temp) * wigner(3) * bessel_plus + exp(m_zi*temp2) * wigner(2) * bessel_minus)
958 end if
959 ! Calculate complex Ay component if generalized bessel OR y -polarized bessel
960 if (pol/=1) then
961 bessel_field(ip, 2) = m_zi * amp * (-exp(m_zi*temp) * wigner(3) * bessel_plus + &
962 exp(m_zi*temp2) * wigner(2) * bessel_minus)
963 end if
964 ! Calculate complex Az component, only iff generalized Bessel
965 if (pol == 0) then
966 bessel_field(ip, 3) = - m_zi * amp * sqrt(m_two) * wigner(1) * loct_bessel(mm, kappa * rho / p_c) * exp(m_zi*temp3)
967 end if
968
969 if (this%envelope(iline)) then
970 k_norm = norm2(k_vector)
971 velocity_time = k_vector * p_c * c_factor * time / k_norm
972 x_prop(:) = pos(:) - velocity_time(:)
973 efield_ip = mxf_envelope_eval(envelope_mxf, x_prop)
974 bessel_field(ip, :) = bessel_field(ip, :) * real(efield_ip, real64)
975 end if
976
977 end do
978
979 end subroutine bessel_beam_function
980
982 subroutine bessel_beam_init(this, nlines, dim)
983 class(bessel_beam_t), intent(out) :: this
984 integer, intent(in) :: nlines
985 integer, intent(in) :: dim
986
987 safe_allocate(this%amp(1: nlines))
988 safe_allocate(this%omega(1:nlines))
989 safe_allocate(this%theta_k(1:nlines))
990 safe_allocate(this%m_order(1:nlines))
991 safe_allocate(this%helicity(1:nlines))
992 safe_allocate(this%shift(1:nlines, 1:dim))
993 safe_allocate(this%envelope(1:nlines))
994 safe_allocate(this%lin_dir(1:nlines))
995 this%amp = m_zero
996 this%omega = m_zero
997 this%theta_k = m_zero
998 this%m_order = m_zero
999 this%helicity = m_zero
1000 this%shift = m_zero
1001 this%lin_dir = m_zero
1002 this%envelope = .false.
1003
1004 end subroutine bessel_beam_init
1005
1007 subroutine bessel_beam_finalize(this)
1008 type(bessel_beam_t), intent(inout) :: this
1009
1010 safe_deallocate_a(this%amp)
1011 safe_deallocate_a(this%omega)
1012 safe_deallocate_a(this%theta_k)
1013 safe_deallocate_a(this%m_order)
1014 safe_deallocate_a(this%helicity)
1015 safe_deallocate_a(this%shift)
1016 safe_deallocate_a(this%lin_dir)
1017 safe_deallocate_a(this%envelope)
1018
1019 end subroutine bessel_beam_finalize
1020
1021end module external_waves_oct_m
1022
1023!! Local Variables:
1024!! mode: f90
1025!! coding: utf-8
1026!! End:
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
subroutine, public accel_release_buffer(this)
Definition: accel.F90:1248
pure logical function, public accel_is_enabled()
Definition: accel.F90:401
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public zderivatives_curl(der, ff, op_ff, ghost_update, set_bc)
apply the curl operator to a vector of mesh functions
subroutine, public dderivatives_curl(der, ff, op_ff, ghost_update, set_bc)
apply the curl operator to a vector of mesh functions
subroutine, public load_external_waves(partners, namespace)
subroutine bessel_beam_function(this, iline, shift, mesh, n_points, time, k_vector, c_factor, envelope_mxf, bessel_field)
. Evaluation of the Bessel beam expression
subroutine, public bessel_source_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
Calculation of Bessel beam from parsed formula.
subroutine external_waves_update_quantity(this, label)
subroutine pw_mx_function_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
Evaluate expression for plane wave that uses predefeined Maxwell function.
subroutine bessel_beam_init(this, nlines, dim)
. Initialization of Bessel beam arrays
subroutine external_waves_copy_quantities_to_interaction(partner, interaction)
subroutine, public external_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
Calculation of external waves from parsed formula.
class(external_waves_t) function, pointer external_waves_constructor(namespace)
subroutine, public external_waves_end(external_waves)
subroutine pw_parsed_evaluation(external_waves, wn, time, mesh, n_points, pw_field)
Evaluate expression for plane wave parsing the provided formula.
subroutine plane_waves_eval(external_waves, time, mesh, type_of_field, out_field_total, der)
Calculation of plane waves from parsed formula.
subroutine get_pw_b_field(external_waves, mesh, pwidx, e_field, b_field)
Calculation of magnetic field for a plane wave.
subroutine external_waves_init_interaction_as_partner(partner, interaction)
subroutine, public external_waves_init(external_waves, namespace)
Here, plane wave is evaluated from analytical formulae on grid.
subroutine bessel_beam_finalize(this)
. Finalize Bessel beam arrays
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:224
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
This module implements the index, used for the mesh points.
Definition: index.F90:122
integer, parameter, public mxll_vec_pot_to_matter
integer, parameter, public mxll_b_field_to_matter
integer, parameter, public mxll_e_field_to_matter
This module defines the abstract interaction_t class, and some auxiliary classes for interactions.
This module defines classes and functions for interaction partners.
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
integer, parameter, public e_field_electric
Definition: lasers.F90:177
integer, parameter, public e_field_vector_potential
Definition: lasers.F90:177
complex(real64) function mxf_envelope_eval(f, x)
Evaluation of envelope itself.
subroutine, public mxf_read(f, namespace, function_name, ierr)
This function initializes "f" from the MXFunctions block.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
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.
Definition: mesh.F90:380
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
character(len=512), private msg
Definition: messages.F90:165
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_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
Some general things and nomenclature:
Definition: par_vec.F90:171
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 defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
Definition: quantity.F90:138
subroutine, public conv_to_c_string(str)
converts to c string
Definition: string.F90:252
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
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 representing derivatives
abstract class for general interaction partners
surrogate interaction class to avoid circular dependencies between modules.
Lorenz force between a systems of particles and an electromagnetic field.
Describes mesh distribution to nodes.
Definition: mesh.F90:186
class to transfer a Maxwell B field to a matter system
class to transfer a Maxwell field to a medium
class to transfer a Maxwell vector potential to a medium
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...
Definition: quantity.F90:171
int true(void)