Octopus
maxwell.F90
Go to the documentation of this file.
1!! Copyright (C) 2019-2020 Franco Bonafe, Heiko Appel, Rene Jestaedt
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module maxwell_oct_m
22 use accel_oct_m
23 use batch_oct_m
28 use debug_oct_m
35 use global_oct_m
36 use grid_oct_m
39 use index_oct_m
44 use io_oct_m
49 use loct_oct_m
51 use math_oct_m
53 use mesh_oct_m
56 use mpi_oct_m
62 use output_oct_m
65 use parser_oct_m
77 use sort_oct_m
78 use space_oct_m
79 use system_oct_m
84 use types_oct_m
85 use unit_oct_m
87
88
89 implicit none
90
91 private
92 public :: &
93 maxwell_t, &
95
96 integer, parameter, public :: &
97 MULTIGRID_MX_TO_MA_EQUAL = 1, &
99
102 type, extends(system_t) :: maxwell_t
103 type(space_t) :: space
104 type(states_mxll_t) :: st
105 type(hamiltonian_mxll_t) :: hm
106 type(grid_t) :: gr
107 type(output_t) :: outp
108 type(multicomm_t) :: mc
109
110 type(mesh_interpolation_t) :: mesh_interpolate
111
112 type(propagator_mxll_t) :: tr_mxll
113 type(td_write_t) :: write_handler
114 type(c_ptr) :: output_handle
115
116 complex(real64), allocatable :: ff_rs_inhom_t1(:,:), ff_rs_inhom_t2(:,:)
117 complex(real64), allocatable :: rs_state_init(:,:)
118 type(time_interpolation_t), pointer :: current_interpolation
119 real(real64) :: bc_bounds(2, 3), dt_bounds(2, 3)
120 integer :: energy_update_iter
121 type(restart_t) :: restart_dump
122 type(restart_t) :: restart
123
124 type(lattice_vectors_t) :: latt
125
126 type(helmholtz_decomposition_t) :: helmholtz
127
128 logical :: write_previous_state = .false.
129
130 contains
131 procedure :: init_interaction => maxwell_init_interaction
132 procedure :: init_parallelization => maxwell_init_parallelization
133 procedure :: initialize => maxwell_initialize
134 procedure :: do_algorithmic_operation => maxwell_do_algorithmic_operation
135 procedure :: is_tolerance_reached => maxwell_is_tolerance_reached
136 procedure :: init_interaction_as_partner => maxwell_init_interaction_as_partner
137 procedure :: copy_quantities_to_interaction => maxwell_copy_quantities_to_interaction
138 procedure :: output_start => maxwell_output_start
139 procedure :: output_write => maxwell_output_write
140 procedure :: output_finish => maxwell_output_finish
141 procedure :: update_interactions_start => maxwell_update_interactions_start
142 procedure :: update_interactions_finish => maxwell_update_interactions_finish
143 procedure :: restart_write_data => maxwell_restart_write_data
144 procedure :: restart_read_data => maxwell_restart_read_data
145 procedure :: update_kinetic_energy => maxwell_update_kinetic_energy
146 procedure :: get_current => maxwell_get_current
147 final :: maxwell_finalize
148 end type maxwell_t
149
150 interface maxwell_t
151 procedure maxwell_constructor
152 end interface maxwell_t
153
154contains
155
156 ! ---------------------------------------------------------
157 function maxwell_constructor(namespace) result(sys)
158 class(maxwell_t), pointer :: sys
159 type(namespace_t), intent(in) :: namespace
160
161 push_sub(maxwell_constructor)
162
163 allocate(sys)
164
165 call maxwell_init(sys, namespace)
166
167 pop_sub(maxwell_constructor)
168 end function maxwell_constructor
169
170
171 ! ---------------------------------------------------------
172 subroutine maxwell_init(this, namespace)
173 class(maxwell_t), intent(inout) :: this
174 type(namespace_t), intent(in) :: namespace
175
176
177 push_sub(maxwell_init)
178
179 call profiling_in("MAXWELL_INIT")
180
181 this%namespace = namespace
182
183 call messages_obsolete_variable(this%namespace, 'SystemName')
184 call messages_print_with_emphasis(msg='Maxwell Simulation Start', namespace=this%namespace)
185 this%space = space_t(this%namespace)
186 if (this%space%is_periodic()) then
187 call messages_not_implemented('Maxwell for periodic systems', namespace=namespace)
188 end if
190 call grid_init_stage_1(this%gr, this%namespace, this%space)
191 call states_mxll_init(this%st, this%namespace, this%space)
192 this%latt = lattice_vectors_t(this%namespace, this%space)
193
194 call this%quantities%add(quantity_t("E field", updated_on_demand = .false.))
195 call this%quantities%add(quantity_t("B field", updated_on_demand = .false.))
196 call this%quantities%add(quantity_t("vector potential", updated_on_demand = .false.))
198 call mesh_interpolation_init(this%mesh_interpolate, this%gr)
200 this%supported_interactions = [linear_medium_to_em_field, current_to_mxll_field]
202
203 call profiling_out("MAXWELL_INIT")
204
206 end subroutine maxwell_init
208 ! ---------------------------------------------------------
209 subroutine maxwell_init_interaction(this, interaction)
210 class(maxwell_t), target, intent(inout) :: this
211 class(interaction_t), intent(inout) :: interaction
213 integer :: depth
216
217 select type (interaction)
219 call interaction%init(this%gr)
221 call interaction%init(this%gr, this%st%dim)
222 call interaction%init_space_latt(this%space, this%latt)
223
224 ! set interpolation depth for interaction
225 ! interpolation depth depends on the propagator
226 select type (prop => this%algo)
228 depth = 2
230 depth = 2
232 depth = 2
234 depth = 4
235 class default
236 message(1) = "The chosen propagator does not yet support interaction interpolation"
237 call messages_fatal(1, namespace=this%namespace)
238 end select
239 call interaction%init_interpolation(depth, interaction%label, cmplx=.true.)
241 this%hm%current_density_from_medium = .true.
242 class default
243 message(1) = "Trying to initialize an unsupported interaction by Maxwell."
244 call messages_fatal(1, namespace=this%namespace)
245 end select
246
248 end subroutine maxwell_init_interaction
249
250 ! ---------------------------------------------------------
251 subroutine maxwell_init_parallelization(this, grp)
252 class(maxwell_t), intent(inout) :: this
253 type(mpi_grp_t), intent(in) :: grp
254
255 integer(int64) :: index_range(4)
256 integer :: ierr, ip, pos_index, rankmin
257 real(real64) :: dmin
258
260
261 call system_init_parallelization(this, grp)
262
263 ! store the ranges for these two indices (serves as initial guess
264 ! for parallelization strategy)
265 index_range(1) = this%gr%np_global ! Number of points in mesh
266 index_range(2) = this%st%nst ! Number of states
267 index_range(3) = 1 ! Number of k-points
268 index_range(4) = 100000 ! Some large number
269
270 ! create index and domain communicators
271 call multicomm_init(this%mc, this%namespace, mpi_world, calc_mode_par, mpi_world%size, &
272 index_range, (/ 5000, 1, 1, 1 /))
273
274 call grid_init_stage_2(this%gr, this%namespace, this%space, this%mc)
275 call output_mxll_init(this%outp, this%namespace, this%space)
276 call hamiltonian_mxll_init(this%hm, this%namespace, this%gr, this%st)
277
278 this%st%energy_rate = m_zero
279 this%st%delta_energy = m_zero
280 this%st%energy_via_flux_calc = m_zero
281 this%st%trans_energy_rate = m_zero
282 this%st%trans_delta_energy = m_zero
283 this%st%trans_energy_via_flux_calc = m_zero
284 this%st%plane_waves_energy_rate = m_zero
285 this%st%plane_waves_delta_energy = m_zero
286 this%st%plane_waves_energy_via_flux_calc = m_zero
287
288 safe_allocate(this%rs_state_init(1:this%gr%np_part, 1:this%st%dim))
289 this%rs_state_init(:,:) = m_z0
290
291 this%energy_update_iter = 1
292
293 call poisson_init(this%st%poisson, this%namespace, this%space, this%gr%der, this%mc, this%gr%stencil)
294
295 call propagator_mxll_init(this%gr, this%namespace, this%st, this%hm, this%tr_mxll)
296 call states_mxll_allocate(this%st, this%gr)
297 call external_current_init(this%st, this%space, this%namespace, this%gr)
298 this%hm%propagation_apply = .true.
299
300 if (parse_is_defined(this%namespace, 'MaxwellIncidentWaves') .and. (this%tr_mxll%bc_plane_waves)) then
301 this%st%rs_state_plane_waves(:,:) = m_z0
302 end if
303
304 ! set map for selected points
305 do ip = 1, this%st%selected_points_number
306 pos_index = mesh_nearest_point(this%gr, this%st%selected_points_coordinate(:,ip), dmin, rankmin)
307 if (this%gr%mpi_grp%rank == rankmin) then
308 this%st%selected_points_map(ip) = pos_index
309 else
310 this%st%selected_points_map(ip) = -1
311 end if
312 end do
313 if (accel_is_enabled()) then
314 call accel_create_buffer(this%st%buff_selected_points_map, accel_mem_read_only, type_integer, &
315 this%st%selected_points_number)
316 call accel_write_buffer(this%st%buff_selected_points_map, this%st%selected_points_number, &
317 this%st%selected_points_map)
318 end if
319
320 this%hm%plane_waves_apply = .true.
321 this%hm%spatial_constant_apply = .true.
322
323 call bc_mxll_init(this%hm%bc, this%namespace, this%space, this%gr, this%st)
324 this%bc_bounds(:,1:3) = this%hm%bc%bc_bounds(:,1:3)
325 call inner_and_outer_points_mapping(this%gr, this%st, this%bc_bounds)
326 this%dt_bounds(2, 1:3) = this%bc_bounds(1, 1:3)
327 this%dt_bounds(1, 1:3) = this%bc_bounds(1, 1:3) - this%gr%der%order * this%gr%spacing(1:3)
328 call surface_grid_points_mapping(this%gr, this%st, this%dt_bounds)
329
330 call restart_init(this%restart, this%namespace, restart_td, restart_type_load, this%mc, ierr, mesh=this%gr)
331 call restart_init(this%restart_dump, this%namespace, restart_td, restart_type_dump, this%mc, ierr, mesh=this%gr)
332
333 ! initialize batches
334 call zbatch_init(this%st%rs_stateb, 1, 1, this%st%dim, this%gr%np_part)
335 if (this%st%pack_states) call this%st%rs_stateb%do_pack()
336 call this%st%rs_stateb%copy_to(this%st%rs_state_prevb)
337 call this%st%rs_stateb%copy_to(this%st%inhomogeneousb)
338 if (this%tr_mxll%bc_plane_waves) then
339 call this%st%rs_stateb%copy_to(this%st%rs_state_plane_wavesb)
340 end if
341
342 ! the order should depend on the propagator
343 !this%current_interpolation => time_interpolation_t(this%gr%np, this%st%dim, 2, .true., "current")
344 ! Initialize Helmholtz decomposition
345 call this%helmholtz%init(this%namespace, this%gr, this%mc, this%space)
346
348 end subroutine maxwell_init_parallelization
349
350 ! ---------------------------------------------------------
351 subroutine maxwell_initialize(this)
352 class(maxwell_t), intent(inout) :: this
353
354 real(real64) :: courant
355
356
357 push_sub(maxwell_initialize)
358
359 call profiling_in("MAXWELL_INIT_CONDITIONS")
360
361 select type (algo => this%algo)
362 class is (propagator_t)
363
364 ! The reference Courant–Friedrichs–Lewy condition for stability here
365 ! is calculated with respect to R. Courant, K. Friedrichs, and H. Lewy
366 ! On the Partial Difference Equations of
367 ! Mathematical Physics, IBM Journal of Research and Development, vol. 11
368 ! pp. 215-234, Mar. 1967.
369 courant = m_one/(p_c * sqrt(m_one/this%gr%spacing(1)**2 + m_one/this%gr%spacing(2)**2 + &
370 m_one/this%gr%spacing(3)**2))
371
372 if (algo%dt > (m_one + m_fourth) * courant) then
373 write(message(1),'(a)') 'The timestep is too large compared to the Courant-Friedrichs-Lewy'
374 write(message(2),'(a)') 'stability condition. Time propagation might not be stable.'
375 call messages_warning(2, namespace=this%namespace)
376 end if
377
378 if (parse_is_defined(this%namespace, 'UserDefinedInitialMaxwellStates')) then
379 call states_mxll_read_user_def(this%namespace, this%space, this%gr, this%gr%der, this%st, this%hm%bc,&
380 this%rs_state_init)
381 call messages_print_with_emphasis(msg="Setting initial EM field inside box", namespace=this%namespace)
382 ! TODO: add consistency check that initial state fulfills Gauss laws
383 this%st%rs_state(:,:) = this%st%rs_state + this%rs_state_init
384 if (this%tr_mxll%bc_plane_waves) then
385 this%st%rs_state_plane_waves(:,:) = this%rs_state_init
386 end if
387 end if
388
389 ! initialize the spatial constant field according to the conditions set in the
390 ! UserDefinedConstantSpatialMaxwellField block
391 if (this%tr_mxll%bc_constant) then
392 call spatial_constant_calculation(this%tr_mxll%bc_constant, this%st, this%gr, this%hm, m_zero, &
393 algo%dt, m_zero, this%st%rs_state, set_initial_state = .true.)
394 ! for mesh parallelization, this needs communication!
395 this%st%rs_state_const(:) = this%st%rs_state(mesh_global_index_from_coords(this%gr, [0,0,0]),:)
396 end if
397
398 if (parse_is_defined(this%namespace, 'UserDefinedInitialMaxwellStates')) then
399 safe_deallocate_a(this%rs_state_init)
400 end if
401
402 call hamiltonian_mxll_update(this%hm, time = m_zero)
403
404 ! calculate Maxwell energy
405 call energy_mxll_calc(this%gr, this%st, this%hm, this%hm%energy, this%st%rs_state, &
406 this%st%rs_state_plane_waves)
407
408 ! Get RS states values for selected points
409 call get_rs_state_at_point(this%st%selected_points_rs_state(:,:), this%st%rs_state, &
410 this%st%selected_points_coordinate(:,:), this%st, this%gr)
411
412 call mxll_set_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
413 call batch_set_zero(this%st%inhomogeneousb)
414 if (this%tr_mxll%bc_plane_waves) then
415 call mxll_set_batch(this%st%rs_state_plane_wavesb, this%st%rs_state_plane_waves, this%gr%np, this%st%dim)
416 end if
417 end select
418
419 call profiling_out("MAXWELL_INIT_CONDITIONS")
420
421 pop_sub(maxwell_initialize)
422 end subroutine maxwell_initialize
423
424 ! ---------------------------------------------------------
425 logical function maxwell_do_algorithmic_operation(this, operation, updated_quantities) result(done)
426 class(maxwell_t), intent(inout) :: this
427 class(algorithmic_operation_t), intent(in) :: operation
428 character(len=:), allocatable, intent(out) :: updated_quantities(:)
429
430 complex(real64), allocatable :: charge_density_ext(:)
431 real(real64) :: current_time
432
434 call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
435
436 current_time = this%iteration%value()
437
438 select type (algo => this%algo)
439 class is (propagator_t)
440
441 done = .true.
442 select case (operation%id)
444 ! For the moment we do nothing
445
446 case (expmid_start)
447 safe_allocate(this%ff_rs_inhom_t1(1:this%gr%np_part, 1:this%hm%dim))
448 safe_allocate(this%ff_rs_inhom_t2(1:this%gr%np_part, 1:this%hm%dim))
449
450 case (expmid_finish)
451 safe_deallocate_a(this%ff_rs_inhom_t1)
452 safe_deallocate_a(this%ff_rs_inhom_t2)
453
454 case (expmid_extrapolate)
455 if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
456 call this%get_current(current_time, this%st%rs_current_density_t1)
457 call this%get_current(current_time+algo%dt, this%st%rs_current_density_t2)
458 end if
459
460 case (expmid_propagate)
461 ! Propagation
462
463 !We first compute three external charge and current densities and we convert them as RS vectors
464 safe_allocate(charge_density_ext(1:this%gr%np))
465
466 !No charge density at the moment
467 charge_density_ext = m_z0
468
469 call transform_rs_densities(this%hm, this%gr, charge_density_ext, &
470 this%st%rs_current_density_t1, this%ff_rs_inhom_t1, rs_trans_forward)
471 call transform_rs_densities(this%hm, this%gr, charge_density_ext, &
472 this%st%rs_current_density_t2, this%ff_rs_inhom_t2, rs_trans_forward)
473
474 safe_deallocate_a(charge_density_ext)
475
476 ! Propagation dt with H_maxwell
477 call mxll_propagation_step(this%hm, this%namespace, this%gr, this%space, this%st, this%tr_mxll,&
478 this%st%rs_stateb, this%ff_rs_inhom_t1, this%ff_rs_inhom_t2, current_time, algo%dt)
479
480 updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
481
482 case (leapfrog_start)
483 if (any(this%hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
484 call bc_mxll_initialize_pml_simple(this%hm%bc, this%space, this%gr, this%hm%c_factor, algo%dt)
485 end if
486
487 case (leapfrog_finish)
488
489 case (leapfrog_propagate)
490 if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
491 call this%get_current(current_time, this%st%rs_current_density_t1)
492 call mxll_set_batch(this%st%inhomogeneousb, this%st%rs_current_density_t1, this%gr%np, this%st%dim)
493 call batch_scal(this%gr%np, -m_one, this%st%inhomogeneousb)
494 end if
495
496 call mxll_propagate_leapfrog(this%hm, this%namespace, this%gr, this%space, this%st, this%tr_mxll, &
497 current_time, algo%dt, this%iteration%counter())
498
499 updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
500
501 case (exp_gauss1_start)
502 safe_allocate(this%ff_rs_inhom_t1(1:this%gr%np_part, 1:this%hm%dim))
503 if (any(this%hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
504 call bc_mxll_initialize_pml_simple(this%hm%bc, this%space, this%gr, this%hm%c_factor, algo%dt)
505 end if
506
507 case (exp_gauss1_finish)
508 safe_deallocate_a(this%ff_rs_inhom_t1)
509
511 if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
512 call this%get_current(current_time+algo%dt*m_half, this%st%rs_current_density_t1)
513 end if
514
516 ! the propagation step is
517 ! F_{n+1} = F_n + dt * phi_1(-i*dt*H) [-i*H F_n - J_1]
518 ! where J_1 = J(t_n + dt/2)
519 call mxll_propagate_expgauss1(this%hm, this%namespace, this%gr, this%space, this%st, this%tr_mxll, &
520 current_time, algo%dt)
521
522 updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
523
524 case (exp_gauss2_start)
525 safe_allocate(this%ff_rs_inhom_t1(1:this%gr%np_part, 1:this%hm%dim))
526 safe_allocate(this%ff_rs_inhom_t2(1:this%gr%np_part, 1:this%hm%dim))
527 if (any(this%hm%bc%bc_ab_type == option__maxwellabsorbingboundaries__cpml)) then
528 call bc_mxll_initialize_pml_simple(this%hm%bc, this%space, this%gr, this%hm%c_factor, algo%dt)
529 end if
530
531 case (exp_gauss2_finish)
532 safe_deallocate_a(this%ff_rs_inhom_t1)
533 safe_deallocate_a(this%ff_rs_inhom_t2)
534
536 if (this%hm%current_density_ext_flag .or. this%hm%current_density_from_medium) then
537 call this%get_current(current_time+algo%dt*(m_half-sqrt(m_three)/6.0_real64), this%st%rs_current_density_t1)
538 call this%get_current(current_time+algo%dt*(m_half+sqrt(m_three)/6.0_real64), this%st%rs_current_density_t2)
539 end if
540
542 ! the propagation step is
543 ! F_{n+1} = F_n + dt * phi_1(-i*dt*H) [-i*H F_n - a_1 J_1 - a_2 J_2]
544 ! + dt * phi_2(-i*dt*H) [-b_1 J_1 - b_2 J_2]
545 ! where J_1 = J(t_n + (1/2 - sqrt(3)/6)*dt),
546 ! J_2 = J(t_n + (1/2 + sqrt(3)/6)*dt),
547 ! a_1 = 1/2*(1+sqrt(3)),
548 ! a_2 = 1/2*(1-sqrt(3)),
549 ! b_1 = -sqrt(3),
550 ! b_2 = sqrt(3)
551 call mxll_propagate_expgauss2(this%hm, this%namespace, this%gr, this%space, this%st, this%tr_mxll, &
552 current_time, algo%dt)
553
554 updated_quantities = [character(16) :: "E field", "B field", "vector potential"]
555
556 case (iteration_done)
558 done = .false.
559
560 case default
561 done = .false.
562 end select
563
564 end select
565
566 call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
569
570 ! ---------------------------------------------------------
571 logical function maxwell_is_tolerance_reached(this, tol) result(converged)
572 class(maxwell_t), intent(in) :: this
573 real(real64), intent(in) :: tol
574
576
577 converged = .false.
578
581
582 ! ---------------------------------------------------------
583 subroutine maxwell_init_interaction_as_partner(partner, interaction)
584 class(maxwell_t), intent(in) :: partner
585 class(interaction_surrogate_t), intent(inout) :: interaction
586
588
589 select type (interaction)
590 type is (lorentz_force_t)
591 ! Nothing to be initialized for the Lorentz force.
593 call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
595 call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
597 call interaction%init_from_partner(partner%gr, partner%space, partner%namespace)
598 class default
599 message(1) = "Unsupported interaction."
600 call messages_fatal(1, namespace=partner%namespace)
601 end select
602
605
606 ! ---------------------------------------------------------
607 subroutine maxwell_copy_quantities_to_interaction(partner, interaction)
608 class(maxwell_t), intent(inout) :: partner
609 class(interaction_surrogate_t), intent(inout) :: interaction
610
611 integer :: ip
612 complex(real64) :: interpolated_value(3)
613 real(real64), allocatable :: b_field(:,:), vec_pot(:,:)
614
616 call profiling_in(trim(partner%namespace%get())//":"//"COPY_QUANTITY_INTER")
617
618 select type (interaction)
619 type is (lorentz_force_t)
620 call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
621
622 do ip = 1, interaction%system_np
623 call mesh_interpolation_evaluate(partner%mesh_interpolate, partner%st%rs_state(:,1), &
624 interaction%system_pos(:, ip), interpolated_value(1))
625 call mesh_interpolation_evaluate(partner%mesh_interpolate, partner%st%rs_state(:,2), &
626 interaction%system_pos(:, ip), interpolated_value(2))
627 call mesh_interpolation_evaluate(partner%mesh_interpolate, partner%st%rs_state(:,3), &
628 interaction%system_pos(:, ip), interpolated_value(3))
629 call get_electric_field_vector(interpolated_value, interaction%partner_e_field(:, ip))
630 call get_magnetic_field_vector(interpolated_value, 1, interaction%partner_b_field(:, ip))
631 end do
632
634 call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
635 select case (interaction%type)
636
637 case (mxll_field_total)
638 call get_electric_field_state(partner%st%rs_state, partner%gr, interaction%partner_field)
639
640 case (mxll_field_trans)
641 call get_transverse_rs_state(partner%helmholtz, partner%st, partner%namespace)
642 call get_electric_field_state(partner%st%rs_state_trans, partner%gr, interaction%partner_field)
643
644 case (mxll_field_long)
645 call partner%helmholtz%get_long_field(partner%namespace, partner%st%rs_state_long, total_field=partner%st%rs_state)
646 call get_electric_field_state(partner%st%rs_state_long, partner%gr, interaction%partner_field)
647
648 case default
649 message(1) = "Unknown type of field requested by interaction."
650 call messages_fatal(1, namespace=partner%namespace)
651 end select
652 call interaction%do_mapping()
653
655 call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
656 safe_allocate(b_field(1:partner%gr%np_part, 1:partner%gr%box%dim))
657 safe_allocate(vec_pot(1:partner%gr%np_part, 1:partner%gr%box%dim))
658 ! magnetic field is always transverse
659 call get_magnetic_field_state(partner%st%rs_state, partner%gr, partner%st%rs_sign, b_field, &
660 partner%st%mu(1:partner%gr%np), partner%gr%np)
661 ! vector potential stored in partner_field
662 call partner%helmholtz%get_vector_potential(partner%namespace, vec_pot, trans_field=b_field)
663 ! in the convention used by the electronic system, the -1/c factor is included in the vector potential
664 ! but we do not divide by it, because the Maxwell code is in SI units, so the vector potential units
665 ! are suitable for a (p + q.A) minimal coupling interaction
666 interaction%partner_field(1:partner%gr%np,1:partner%gr%box%dim) = &
667 vec_pot(1:partner%gr%np,1:partner%gr%box%dim)
668 safe_deallocate_a(b_field)
669 safe_deallocate_a(vec_pot)
670 call interaction%do_mapping()
671
673 call mxll_get_batch(partner%st%rs_stateb, partner%st%rs_state, partner%gr%np, partner%st%dim)
674 call get_magnetic_field_state(partner%st%rs_state, partner%gr, partner%st%rs_sign, &
675 interaction%partner_field, partner%st%mu(1:partner%gr%np), partner%gr%np)
676 call interaction%do_mapping()
677
678 class default
679 message(1) = "Unsupported interaction."
680 call messages_fatal(1, namespace=partner%namespace)
681 end select
682
683 call profiling_out(trim(partner%namespace%get())//":"//"COPY_QUANTITY_INTER")
686
687 ! ---------------------------------------------------------
688 subroutine maxwell_output_start(this)
689 class(maxwell_t), intent(inout) :: this
690
691 real(real64) :: itr_value
692
693 push_sub(maxwell_output_start)
694 call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_START")
695
696 select type (algo => this%algo)
697 class is (propagator_t)
698 call mxll_get_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
699
700 call td_write_mxll_init(this%write_handler, this%namespace, this%iteration%counter(), algo%dt)
701 if (this%st%fromScratch) then
702 call td_write_mxll_iter(this%write_handler, this%space, this%gr, this%st, this%hm, this%helmholtz, algo%dt, &
703 this%iteration%counter(), this%namespace)
704 itr_value = this%iteration%value()
705 call td_write_mxll_free_data(this%write_handler, this%namespace, this%space, this%gr, this%st, this%hm, this%helmholtz, &
706 this%outp, this%iteration%counter(), itr_value)
707 end if
708
709 end select
710
711 call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_START")
712 pop_sub(maxwell_output_start)
713 end subroutine maxwell_output_start
714
715 ! ---------------------------------------------------------
716 subroutine maxwell_output_write(this)
717 class(maxwell_t), intent(inout) :: this
718
719 logical :: stopping, reached_output_interval
720 real(real64) :: itr_value
721 integer :: iout
722
723 push_sub(maxwell_output_write)
724 call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
725
726 select type (algo => this%algo)
727 class is (propagator_t)
728 stopping = clean_stop(this%mc%master_comm)
729
730 call td_write_mxll_iter(this%write_handler, this%space, this%gr, this%st, this%hm, this%helmholtz, algo%dt, &
731 this%iteration%counter(), this%namespace)
732
733 reached_output_interval = .false.
734 do iout = 1, out_maxwell_max
735 if (this%outp%output_interval(iout) > 0) then
736 if (mod(this%iteration%counter(), this%outp%output_interval(iout)) == 0) then
737 reached_output_interval = .true.
738 exit
739 end if
740 end if
741 end do
742
743 if (reached_output_interval .or. stopping) then
744 call mxll_get_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
745
746 itr_value = this%iteration%value()
747 call td_write_mxll_free_data(this%write_handler, this%namespace, this%space, this%gr, this%st, this%hm, this%helmholtz, &
748 this%outp, this%iteration%counter(), itr_value)
749 end if
750
751 end select
752
753 call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
754 pop_sub(maxwell_output_write)
755 end subroutine maxwell_output_write
756
757 ! ---------------------------------------------------------
758 subroutine maxwell_output_finish(this)
759 class(maxwell_t), intent(inout) :: this
760
761
762 push_sub(maxwell_output_finish)
763
764 call profiling_in("MAXWELL_OUTPUT_FINISH")
765
766 select type (algo => this%algo)
767 class is (propagator_t)
768 call td_write_mxll_end(this%write_handler)
769 end select
770
771 call profiling_out("MAXWELL_OUTPUT_FINISH")
772
773 pop_sub(maxwell_output_finish)
774 end subroutine maxwell_output_finish
775
776 ! ---------------------------------------------------------
777 subroutine maxwell_update_interactions_start(this)
778 class(maxwell_t), intent(inout) :: this
779
780 type(interaction_iterator_t) :: iter
781 integer :: int_counter
782
784
785 int_counter = 0
786 call iter%start(this%interactions)
787 do while (iter%has_next())
788 select type (interaction => iter%get_next())
790 int_counter = int_counter + 1
791 end select
792 end do
793
794 if (int_counter /= 0 .and. .not. allocated(this%hm%medium_boxes)) then
795 safe_allocate(this%hm%medium_boxes(1:int_counter))
796 this%hm%calc_medium_box = .true.
797 end if
798
801
802 ! ---------------------------------------------------------
804 class(maxwell_t), intent(inout) :: this
805
806 type(interaction_iterator_t) :: iter
807 integer :: iint
808
810
811 iint = 0
812
813 call iter%start(this%interactions)
814 do while (iter%has_next())
815 select type (interaction => iter%get_next())
817 if (allocated(this%hm%medium_boxes) .and. .not. this%hm%medium_boxes_initialized) then
818 iint = iint + 1
819 this%hm%medium_boxes(iint) = interaction%medium_box
820 end if
821 end select
822 end do
823
824 if (allocated(this%hm%medium_boxes) .and. .not. this%hm%medium_boxes_initialized) then
825 call set_medium_rs_state(this%st, this%gr, this%hm)
826 this%hm%medium_boxes_initialized = .true.
827 end if
828
829 if (this%hm%medium_boxes_initialized .and. this%hm%operator == faraday_ampere) then
830 message(1) = "A linear medium has been defined in the input file but the Hamiltonian"
831 message(2) = "type you specified is not capable of dealing with the medium."
832 message(3) = "Please use MaxwellHamiltonianOperator = faraday_ampere_medium or simple to enable"
833 message(4) = "the medium propagation."
834 call messages_fatal(4, namespace=this%namespace)
835 end if
836
837 if (.not. this%hm%medium_boxes_initialized .and. this%hm%operator == faraday_ampere_medium) then
838 message(1) = "The variable MaxwellHamiltonianOperator has been defined as faraday_ampere_medium"
839 message(2) = "in the input file but no linear medium has been defined in the system block."
840 message(3) = "Please either use a different option for MaxwellHamiltonianOperator or add"
841 message(4) = "a linear medium to the system block."
842 call messages_fatal(4, namespace=this%namespace)
843 end if
844
847
848 ! ---------------------------------------------------------
849 subroutine maxwell_restart_write_data(this)
850 class(maxwell_t), intent(inout) :: this
852 integer :: ierr, err, zff_dim, id, id1, id2, ip_in, offset, iout
853 logical :: pml_check, write_previous_state
854 complex(real64), allocatable :: zff(:,:)
855 type(interaction_iterator_t) :: iter
856 class(interaction_t), pointer :: interaction
857
859 call profiling_in(trim(this%namespace%get())//":"//"RESTART_WRITE")
860
861 ierr = 0
862
863 if (mpi_grp_is_root(mpi_world)) then
864 do iout = 1, out_maxwell_max
865 if (this%write_handler%out(iout)%write) then
866 call write_iter_flush(this%write_handler%out(iout)%handle)
867 end if
868 end do
869 end if
871 if (.not. restart_skip(this%restart_dump)) then
872 pml_check = any(this%hm%bc%bc_ab_type(1:3) == option__maxwellabsorbingboundaries__cpml)
873
874 message(1) = "Debug: Writing td_maxwell restart."
875 call messages_info(1, namespace=this%namespace, debug_only=.true.)
876
877 if (this%tr_mxll%bc_plane_waves) then
878 zff_dim = 2 * this%st%dim
879 else
880 zff_dim = 1 * this%st%dim
881 end if
882 if (pml_check) then
883 zff_dim = zff_dim + 18
884 end if
885 select type (prop => this%algo)
886 type is (propagator_leapfrog_t)
887 write_previous_state = .true.
888 zff_dim = zff_dim + this%st%dim
889 class default
890 write_previous_state = .false.
891 end select
892 if (pml_check .and. accel_is_enabled()) then
893 call accel_read_buffer(this%hm%bc%pml%buff_conv_plus, &
894 int(this%hm%bc%pml%points_number, int64)*3*3, this%hm%bc%pml%conv_plus)
895 call accel_read_buffer(this%hm%bc%pml%buff_conv_minus, &
896 int(this%hm%bc%pml%points_number, int64)*3*3, this%hm%bc%pml%conv_minus)
897 end if
898
899 call mxll_get_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
900 if (write_previous_state) then
901 call mxll_get_batch(this%st%rs_state_prevb, this%st%rs_state_prev, this%gr%np, this%st%dim)
902 end if
903
904 safe_allocate(zff(1:this%gr%np,1:zff_dim))
905 zff = m_z0
906
907 zff(1:this%gr%np, 1:this%st%dim) = this%st%rs_state(1:this%gr%np, 1:this%st%dim)
908 if (this%tr_mxll%bc_plane_waves) then
909 call mxll_get_batch(this%st%rs_state_plane_wavesb, this%st%rs_state_plane_waves, &
910 this%gr%np, this%st%dim)
911 zff(1:this%gr%np, this%st%dim+1:this%st%dim+this%st%dim) = &
912 this%st%rs_state_plane_waves(1:this%gr%np, 1:this%st%dim)
913 offset = 2*this%st%dim
914 else
915 offset = this%st%dim
916 end if
917 if (pml_check) then
918 id = 0
919 do id1 = 1, 3
920 do id2 = 1, 3
921 id = id + 1
922 do ip_in = 1, this%hm%bc%pml%points_number
923 zff(ip_in, offset+id) = this%hm%bc%pml%conv_plus(ip_in, id1, id2)
924 zff(ip_in, offset+9+id) = this%hm%bc%pml%conv_minus(ip_in, id1, id2)
925 end do
926 end do
927 end do
928 offset = offset + 18
929 end if
930 if (write_previous_state) then
931 zff(1:this%gr%np, offset+1:offset+this%st%dim) = &
932 this%st%rs_state_prev(1:this%gr%np, 1:this%st%dim)
933 end if
934
935 call states_mxll_dump(this%restart_dump, this%st, this%space, this%gr, zff, zff_dim, err, this%iteration%counter())
936 if (err /= 0) ierr = ierr + 1
937
938 if (this%hm%current_density_from_medium) then
939 !call this%current_interpolation%write_restart(this%gr, this%space, this%restart_dump, err)
940 call iter%start(this%interactions)
941 do while (iter%has_next())
942 interaction => iter%get_next()
943 select type (interaction)
944 class is (current_to_mxll_field_t)
945 call interaction%write_restart(this%gr, this%space, this%restart_dump, err)
946 end select
947 end do
948 if (err /= 0) ierr = ierr + 1
949 end if
950
951 message(1) = "Debug: Writing td_maxwell restart done."
952 call messages_info(1, namespace=this%namespace, debug_only=.true.)
953
954 safe_deallocate_a(zff)
955
956 if (ierr /=0) then
957 message(1) = "Unable to write time-dependent Maxwell restart information."
958 call messages_warning(1, namespace=this%namespace)
959 end if
960 end if
961
962 call profiling_out(trim(this%namespace%get())//":"//"RESTART_WRITE")
964 end subroutine maxwell_restart_write_data
965
966 ! ---------------------------------------------------------
967 ! this function returns true if restart data could be read
968 logical function maxwell_restart_read_data(this)
969 class(maxwell_t), intent(inout) :: this
970
971 integer :: ierr, err, zff_dim, id, id1, id2, ip_in, offset
972 logical :: pml_check, read_previous_state
973 complex(real64), allocatable :: zff(:,:)
974 type(interaction_iterator_t) :: iter
975 class(interaction_t), pointer :: interaction
976
978 call profiling_in(trim(this%namespace%get())//":"//"RESTART_READ")
979
980 if (.not. restart_skip(this%restart)) then
981 ierr = 0
982 pml_check = any(this%hm%bc%bc_ab_type(1:3) == option__maxwellabsorbingboundaries__cpml)
983
984 if (restart_skip(this%restart)) then
985 ierr = -1
986 call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
988 return
989 end if
990
991 message(1) = "Debug: Reading td_maxwell restart."
992 call messages_info(1, namespace=this%namespace, debug_only=.true.)
993
994 if (this%tr_mxll%bc_plane_waves) then
995 zff_dim = 2 * this%st%dim
996 else
997 zff_dim = 1 * this%st%dim
998 end if
999 if (pml_check) then
1000 zff_dim = zff_dim + 18
1001 end if
1002 select type (prop => this%algo)
1003 type is (propagator_leapfrog_t)
1004 read_previous_state = .true.
1005 zff_dim = zff_dim + this%st%dim
1006 class default
1007 read_previous_state = .false.
1008 end select
1009
1010 safe_allocate(zff(1:this%gr%np,1:zff_dim))
1011
1012 call states_mxll_load(this%restart, this%st, this%gr, this%namespace, this%space, zff, &
1013 zff_dim, err, label = ": td_maxwell")
1014 this%st%rs_current_density_restart = .true.
1015
1016 this%st%rs_state(1:this%gr%np,1:this%st%dim) = zff(1:this%gr%np, 1:this%st%dim)
1017 if (this%tr_mxll%bc_plane_waves) then
1018 this%st%rs_state_plane_waves(1:this%gr%np,1:this%st%dim) = &
1019 zff(1:this%gr%np,this%st%dim+1:this%st%dim+3)
1020 offset = 2*this%st%dim
1021 else
1022 offset = this%st%dim
1023 end if
1024 if (pml_check) then
1025 id = 0
1026 do id1 = 1, 3
1027 do id2 = 1, 3
1028 id = id+1
1029 do ip_in = 1, this%hm%bc%pml%points_number
1030 this%hm%bc%pml%conv_plus(ip_in,id1,id2) = zff(ip_in, offset+ id)
1031 this%hm%bc%pml%conv_minus(ip_in,id1,id2) = zff(ip_in, offset+9+id)
1032 end do
1033 end do
1034 end do
1035 this%hm%bc%pml%conv_plus_old = this%hm%bc%pml%conv_plus
1036 this%hm%bc%pml%conv_minus_old = this%hm%bc%pml%conv_minus
1037 offset = offset + 18
1038 end if
1039 if (read_previous_state) then
1040 this%st%rs_state_prev(1:this%gr%np, 1:this%st%dim) = &
1041 zff(1:this%gr%np, offset+1:offset+this%st%dim)
1042 end if
1043
1044 if (err /= 0) then
1045 ierr = ierr + 1
1046 end if
1047
1048 message(1) = "Debug: Reading td restart done."
1049 call messages_info(1, namespace=this%namespace, debug_only=.true.)
1050
1051 safe_deallocate_a(zff)
1052
1053 if (pml_check .and. accel_is_enabled()) then
1054 call accel_write_buffer(this%hm%bc%pml%buff_conv_plus, &
1055 int(this%hm%bc%pml%points_number, int64)*3*3, this%hm%bc%pml%conv_plus)
1056 call accel_write_buffer(this%hm%bc%pml%buff_conv_minus, &
1057 int(this%hm%bc%pml%points_number, int64)*3*3, this%hm%bc%pml%conv_minus)
1058 call accel_write_buffer(this%hm%bc%pml%buff_conv_plus_old, &
1059 int(this%hm%bc%pml%points_number, int64)*3*3, this%hm%bc%pml%conv_plus_old)
1060 end if
1062 if (this%hm%current_density_from_medium) then
1063 !call this%current_interpolation%read_restart(this%gr, this%space, this%restart, err)
1064 call iter%start(this%interactions)
1065 do while (iter%has_next())
1066 interaction => iter%get_next()
1067 select type (interaction)
1068 class is (current_to_mxll_field_t)
1069 call interaction%read_restart(this%gr, this%space, this%restart, err)
1070 end select
1071 end do
1072 if (err /= 0) then
1073 ierr = ierr + 1
1074 end if
1075 end if
1076
1077 ! set batches
1078 call mxll_set_batch(this%st%rs_stateb, this%st%rs_state, this%gr%np, this%st%dim)
1079 if (read_previous_state) then
1080 call mxll_set_batch(this%st%rs_state_prevb, this%st%rs_state_prev, this%gr%np, this%st%dim)
1081 end if
1082 call batch_set_zero(this%st%inhomogeneousb)
1083 if (this%tr_mxll%bc_plane_waves) then
1084 call mxll_set_batch(this%st%rs_state_plane_wavesb, this%st%rs_state_plane_waves, this%gr%np, this%st%dim)
1085 end if
1086
1087 this%st%fromScratch = .false.
1089 else
1090 message(1) = "Unable to read time-dependent Maxwell restart information: Starting from scratch"
1091 call messages_warning(1, namespace=this%namespace)
1093 end if
1094
1095 call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
1097 end function maxwell_restart_read_data
1098
1099 subroutine maxwell_update_kinetic_energy(this)
1100 class(maxwell_t), intent(inout) :: this
1101
1103
1104 ! the energy has already been computed at the end of the timestep
1105
1106 ! the energy of the EM wave is computed and stored in energy_mxll%energy;
1107 ! energy_mxll%energy = energy_mxll%e_energy + energy_mxll%b_energy
1108 ! here this%hm%energy is 'energy_mxll'
1109 this%kinetic_energy = this%hm%energy%energy
1110
1112 end subroutine maxwell_update_kinetic_energy
1113
1114 ! ---------------------------------------------------------
1115 subroutine maxwell_exec_end_of_timestep_tasks(this)
1116 class(maxwell_t), intent(inout) :: this
1117
1119 call profiling_in(trim(this%namespace%get())//":"//"END_TIMESTEP")
1120
1121 ! calculate Maxwell energy
1122 call energy_mxll_calc_batch(this%gr, this%st, this%hm, this%hm%energy, this%st%rs_stateb, this%st%rs_state_plane_wavesb)
1123
1124 ! get RS state values for selected points
1125 call get_rs_state_batch_selected_points(this%st%selected_points_rs_state(:,:), this%st%rs_stateb, &
1126 this%st, this%gr)
1127
1128 call profiling_out(trim(this%namespace%get())//":"//"END_TIMESTEP")
1131
1132 ! ---------------------------------------------------------
1133 subroutine maxwell_get_current(this, time, current)
1134 class(maxwell_t), intent(inout) :: this
1135 real(real64), intent(in) :: time
1136 complex(real64), contiguous, intent(inout) :: current(:, :)
1137
1138 type(interaction_iterator_t) :: iter
1139 complex(real64), allocatable :: current_density_ext(:, :)
1140
1141 push_sub(maxwell_get_current)
1142
1143 safe_allocate(current_density_ext(1:this%gr%np, 1:this%st%dim))
1144 current = m_z0
1145 if (this%hm%current_density_from_medium) then
1146 ! interpolate current from interaction
1147 call iter%start(this%interactions)
1148 do while (iter%has_next())
1149 select type (interaction => iter%get_next())
1150 class is (current_to_mxll_field_t)
1151 call interaction%interpolate(time, current_density_ext)
1152 call lalg_axpy(this%gr%np, 3, m_one, current_density_ext * this%hm%current_factor, current)
1153 end select
1154 end do
1155 end if
1156 ! calculation of external RS density
1157 if (this%hm%current_density_ext_flag) then
1158 call get_rs_density_ext(this%st, this%space, this%gr, time, current_density_ext)
1159 call lalg_axpy(this%gr%np, 3, m_one, current_density_ext, current)
1160 end if
1161 safe_deallocate_a(current_density_ext)
1162
1163 pop_sub(maxwell_get_current)
1164 end subroutine maxwell_get_current
1165
1166 ! ---------------------------------------------------------
1167 subroutine maxwell_finalize(this)
1168 type(maxwell_t), intent(inout) :: this
1169
1170
1171 push_sub(maxwell_finalize)
1172
1173 call profiling_in("MAXWELL_FINALIZE")
1174
1175 call system_end(this)
1176
1177 ! free memory
1178 safe_deallocate_a(this%rs_state_init)
1179 call hamiltonian_mxll_end(this%hm)
1180
1181 call multicomm_end(this%mc)
1182
1183 call this%st%rs_stateb%end()
1184 call this%st%rs_state_prevb%end()
1185 call this%st%inhomogeneousb%end()
1186 if (this%tr_mxll%bc_plane_waves) then
1187 call this%st%rs_state_plane_wavesb%end()
1188 end if
1189
1190 call states_mxll_end(this%st)
1191
1192 call grid_end(this%gr)
1193
1194 call restart_end(this%restart)
1195 call restart_end(this%restart_dump)
1196
1197 call poisson_end(this%st%poisson)
1198
1199 call profiling_out("MAXWELL_FINALIZE")
1200
1201 pop_sub(maxwell_finalize)
1202 end subroutine maxwell_finalize
1203
1204end module maxwell_oct_m
1205
1206!! Local Variables:
1207!! mode: f90
1208!! coding: utf-8
1209!! End:
scale a batch by a constant or vector
Definition: batch_ops.F90:162
constant times a vector plus a vector
Definition: lalg_basic.F90:171
pure logical function, public accel_is_enabled()
Definition: accel.F90:400
integer, parameter, public accel_mem_read_only
Definition: accel.F90:183
This module defines the abstract interfact for algorithm factories.
This module implements the basic elements defining algorithms.
Definition: algorithm.F90:141
character(len=algo_label_len), parameter, public iteration_done
Definition: algorithm.F90:172
This module implements batches of mesh functions.
Definition: batch.F90:133
subroutine, public zbatch_init(this, dim, st_start, st_end, np, special, packed)
initialize a TYPE_CMPLX valued batch to given size without providing external memory
Definition: batch.F90:1774
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:116
subroutine, public batch_set_zero(this, np, async)
fill all mesh functions of the batch with zero
Definition: batch_ops.F90:242
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public get_rs_density_ext(st, space, mesh, time, rs_current_density_ext)
subroutine, public external_current_init(st, space, namespace, mesh)
This module implements the field transfer.
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_fourth
Definition: global.F90:197
complex(real64), parameter, public m_z0
Definition: global.F90:198
real(real64), parameter, public m_half
Definition: global.F90:194
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
real(real64), parameter, public m_three
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public grid_init_stage_1(gr, namespace, space, symm, latt, n_sites, site_position)
First stage of the grid initialization.
Definition: grid.F90:194
subroutine, public grid_init_stage_2(gr, namespace, space, mc, qvector)
Second stage of the grid initialization.
Definition: grid.F90:465
subroutine, public grid_end(gr)
finalize a grid object
Definition: grid.F90:493
subroutine, public hamiltonian_mxll_init(hm, namespace, gr, st)
Initializing the Maxwell Hamiltonian.
integer, parameter, public faraday_ampere
subroutine, public hamiltonian_mxll_end(hm)
integer, parameter, public faraday_ampere_medium
subroutine, public hamiltonian_mxll_update(this, time)
Maxwell Hamiltonian update (here only the time is updated, can maybe be added to another routine)
The Helmholtz decomposition is intended to contain "only mathematical" functions and procedures to co...
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 linear_medium_to_em_field
integer, parameter, public lorentz_force
integer, parameter, public mxll_b_field_to_matter
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.
Definition: io.F90:114
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public bc_mxll_init(bc, namespace, space, gr, st)
subroutine, public surface_grid_points_mapping(mesh, st, bounds)
subroutine, public bc_mxll_initialize_pml_simple(bc, space, gr, c_factor, dt)
subroutine, public inner_and_outer_points_mapping(mesh, st, bounds)
subroutine maxwell_update_interactions_start(this)
Definition: maxwell.F90:871
subroutine maxwell_exec_end_of_timestep_tasks(this)
Definition: maxwell.F90:1209
subroutine maxwell_init_interaction(this, interaction)
Definition: maxwell.F90:303
subroutine maxwell_output_start(this)
Definition: maxwell.F90:782
subroutine maxwell_update_interactions_finish(this)
Definition: maxwell.F90:897
subroutine maxwell_init_interaction_as_partner(partner, interaction)
Definition: maxwell.F90:677
subroutine maxwell_output_finish(this)
Definition: maxwell.F90:852
subroutine maxwell_restart_write_data(this)
Definition: maxwell.F90:943
logical function maxwell_do_algorithmic_operation(this, operation, updated_quantities)
Definition: maxwell.F90:519
integer, parameter, public multigrid_mx_to_ma_large
Definition: maxwell.F90:189
subroutine maxwell_init_parallelization(this, grp)
Definition: maxwell.F90:345
subroutine maxwell_get_current(this, time, current)
Definition: maxwell.F90:1227
subroutine, public maxwell_init(this, namespace)
Definition: maxwell.F90:266
subroutine maxwell_initialize(this)
Definition: maxwell.F90:445
subroutine maxwell_output_write(this)
Definition: maxwell.F90:810
logical function maxwell_is_tolerance_reached(this, tol)
Definition: maxwell.F90:665
class(maxwell_t) function, pointer maxwell_constructor(namespace)
Definition: maxwell.F90:251
subroutine maxwell_update_kinetic_energy(this)
Definition: maxwell.F90:1193
logical function maxwell_restart_read_data(this)
Definition: maxwell.F90:1062
subroutine maxwell_finalize(this)
Definition: maxwell.F90:1261
subroutine maxwell_copy_quantities_to_interaction(partner, interaction)
Definition: maxwell.F90:701
subroutine, public mesh_interpolation_init(this, mesh)
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
integer(int64) function, public mesh_global_index_from_coords(mesh, ix)
This function returns the true global index of the point for a given vector of integer coordinates.
Definition: mesh.F90:916
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
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
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
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:430
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
subroutine, public multicomm_end(mc)
Definition: multicomm.F90:889
subroutine, public multicomm_init(mc, namespace, base_grp, mode_para, n_node, index_range, min_range)
create index and domain communicators
Definition: multicomm.F90:264
integer, parameter, public mxll_field_total
integer, parameter, public mxll_field_trans
integer, parameter, public mxll_field_long
this module contains the low-level part of the output system
Definition: output_low.F90:115
subroutine, public output_mxll_init(outp, namespace, space)
this module contains the output system
Definition: output.F90:115
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
subroutine, public poisson_init(this, namespace, space, der, mc, stencil, qtot, label, solver, verbose, force_serial, force_cmplx)
Definition: poisson.F90:233
subroutine, public poisson_end(this)
Definition: poisson.F90:709
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
character(len=algo_label_len), parameter, public exp_gauss1_start
character(len=algo_label_len), parameter, public exp_gauss1_finish
character(len=algo_label_len), parameter, public exp_gauss1_extrapolate
character(len=algo_label_len), parameter, public exp_gauss1_propagate
character(len=algo_label_len), parameter, public exp_gauss2_finish
character(len=algo_label_len), parameter, public exp_gauss2_propagate
character(len=algo_label_len), parameter, public exp_gauss2_extrapolate
character(len=algo_label_len), parameter, public exp_gauss2_start
character(len=algo_label_len), parameter, public expmid_extrapolate
character(len=algo_label_len), parameter, public expmid_finish
character(len=algo_label_len), parameter, public expmid_start
character(len=algo_label_len), parameter, public expmid_propagate
character(len=algo_label_len), parameter, public leapfrog_propagate
character(len=algo_label_len), parameter, public leapfrog_finish
character(len=algo_label_len), parameter, public leapfrog_start
subroutine, public mxll_propagate_expgauss2(hm, namespace, gr, space, st, tr, time, dt)
Exponential propagation scheme with Gauss collocation points, s=2.
subroutine, public energy_mxll_calc_batch(gr, st, hm, energy_mxll, rs_fieldb, rs_field_plane_wavesb)
integer, parameter, public rs_trans_forward
subroutine, public mxll_propagate_leapfrog(hm, namespace, gr, space, st, tr, time, dt, counter)
subroutine, public mxll_propagate_expgauss1(hm, namespace, gr, space, st, tr, time, dt)
Exponential propagation scheme with Gauss collocation points, s=1.
subroutine, public mxll_propagation_step(hm, namespace, gr, space, st, tr, rs_stateb, ff_rs_inhom_t1, ff_rs_inhom_t2, time, dt)
subroutine, public spatial_constant_calculation(constant_calc, st, gr, hm, time, dt, delay, rs_state, set_initial_state)
subroutine, public set_medium_rs_state(st, gr, hm)
subroutine, public transform_rs_densities(hm, mesh, rs_charge_density, rs_current_density, ff_density, sign)
subroutine, public energy_mxll_calc(gr, st, hm, energy_mxll, rs_field, rs_field_plane_waves)
subroutine, public propagator_mxll_init(gr, namespace, st, hm, tr)
This module implements the basic propagator framework.
Definition: propagator.F90:117
character(len=algo_label_len), parameter, public store_current_status
Definition: propagator.F90:169
This module defines the quantity_t class and the IDs for quantities, which can be exposed by a system...
Definition: quantity.F90:138
Implementation details for regridding.
Definition: regridding.F90:170
logical function, public clean_stop(comm)
returns true if a file named stop exists
Definition: restart.F90:276
subroutine, public restart_init(restart, namespace, data_type, type, mc, ierr, mesh, dir, exact)
Initializes a restart object.
Definition: restart.F90:516
integer, parameter, public restart_type_dump
Definition: restart.F90:245
logical pure function, public restart_skip(restart)
Returns true if the restart information should neither be read nor written. This might happen because...
Definition: restart.F90:969
integer, parameter, public restart_td
Definition: restart.F90:200
integer, parameter, public restart_type_load
Definition: restart.F90:245
subroutine, public restart_end(restart)
Definition: restart.F90:722
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
subroutine, public mxll_set_batch(rs_stateb, rs_state, np, dim, offset)
subroutine, public get_rs_state_at_point(rs_state_point, rs_state, pos, st, mesh)
subroutine, public get_electric_field_vector(rs_state_vector, electric_field_vector, ep_element)
subroutine, public get_transverse_rs_state(helmholtz, st, namespace)
subroutine, public mxll_get_batch(rs_stateb, rs_state, np, dim, offset)
subroutine, public get_rs_state_batch_selected_points(rs_state_point, rs_stateb, st, mesh)
subroutine, public states_mxll_end(st)
subroutine, public states_mxll_allocate(st, mesh)
Allocates the Maxwell states defined within a states_mxll_t structure.
subroutine, public states_mxll_init(st, namespace, space)
subroutine, public get_electric_field_state(rs_state, mesh, electric_field, ep_field, np)
subroutine, public get_magnetic_field_state(rs_state, mesh, rs_sign, magnetic_field, mu_field, np)
subroutine, public get_magnetic_field_vector(rs_state_vector, rs_sign, magnetic_field_vector, mu_element)
subroutine, public states_mxll_load(restart, st, mesh, namespace, space, zff, zff_dim, ierr, iter, lowest_missing, label, verbose)
subroutine, public states_mxll_dump(restart, st, space, mesh, zff, zff_dim, ierr, iter, st_start_writing, verbose)
subroutine, public states_mxll_read_user_def(namespace, space, mesh, der, st, bc, user_def_rs_state)
This module implements the abstract system type.
Definition: system.F90:118
subroutine, public system_init_parallelization(this, grp)
Basic functionality: copy the MPI group. This function needs to be implemented by extended types that...
Definition: system.F90:1236
subroutine, public system_end(this)
Definition: system.F90:1145
integer, parameter, public out_maxwell_max
Definition: td_write.F90:289
subroutine, public td_write_mxll_init(writ, namespace, iter, dt)
Definition: td_write.F90:3554
subroutine, public td_write_mxll_end(writ)
Definition: td_write.F90:3742
subroutine, public td_write_mxll_free_data(writ, namespace, space, gr, st, hm, helmholtz, outp, iter, time)
Definition: td_write.F90:4200
subroutine, public td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
Definition: td_write.F90:3758
type(type_t), public type_integer
Definition: types.F90:135
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.
Descriptor of one algorithmic operation.
Definition: algorithm.F90:163
Class to transfer a current to a Maxwell field.
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.
Lorenz force between a systems of particles and an electromagnetic field.
Class describing Maxwell systems.
Definition: maxwell.F90:195
This is defined even when running serial.
Definition: mpi.F90:142
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
Implements the an exponential RK scheme with Gauss collocation points, s=1 see also Hochbruck,...
Implements the an exponential RK scheme with Gauss collocation points, s=2 see also Hochbruck,...
Implements the explicit exponential midpoint propagator (without predictor-corrector)
Implements a propagator for the leap frog algorithm.
Abstract class implementing propagators.
Definition: propagator.F90:138
Systems (system_t) can expose quantities that can be used to calculate interactions with other system...
Definition: quantity.F90:171
Abstract class for systems.
Definition: system.F90:172
int true(void)