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