Octopus
external_densities.F90
Go to the documentation of this file.
1!! Copyright (C) 2019 R. Jestaedt, H. Appel
2!! Copyright (C) 2020 N. Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
23 use debug_oct_m
24 use global_oct_m
26 use math_oct_m
29 use mesh_oct_m
31 use mpi_oct_m
33 use parser_oct_m
35 use space_oct_m
37 use types_oct_m
38 use unit_oct_m
42 use xc_oct_m
43
44 implicit none
45
46 private
47
48 public :: &
52
53 integer, parameter, public :: &
54 EXTERNAL_CURRENT_PARSER = 0, &
56
57contains
58
59 !----------------------------------------------------------
60 subroutine get_rs_density_ext(st, space, mesh, time, rs_current_density_ext)
61 type(states_mxll_t), intent(inout) :: st
62 class(space_t), intent(in) :: space
63 class(mesh_t), intent(in) :: mesh
64 real(real64), intent(in) :: time
65 complex(real64), optional, intent(inout) :: rs_current_density_ext(:,:)
66
67 real(real64), allocatable :: current(:,:,:)
68
69 push_sub(get_rs_density_ext)
70
71 safe_allocate(current(1:mesh%np, 1:space%dim, 1))
72
73 call external_current_calculation(st, space, mesh, time, current(:, :, 1))
74 call build_rs_current_state(current(:, :, 1), mesh, rs_current_density_ext, st%ep, mesh%np)
75
76 safe_deallocate_a(current)
77
78 pop_sub(get_rs_density_ext)
79 end subroutine get_rs_density_ext
80
81
82 !----------------------------------------------------------
83 subroutine external_current_init(st, space, namespace, mesh)
84 type(states_mxll_t), intent(inout) :: st
85 class(space_t), intent(in) :: space
86 class(mesh_t), intent(in) :: mesh
87 type(namespace_t), intent(in) :: namespace
88
89 type(block_t) :: blk
90 integer :: ip, il, nlines, ncols, idir, ierr
91 real(real64) :: j_vector(space%dim), dummy(space%dim), xx(1:space%dim), rr, omega
92 character(len=1024) :: tdf_expression, phase_expression
93
94
95 push_sub(external_current_init)
96
97 call profiling_in('EXTERNAL_CURRENT_INIT')
98
99 !%Variable UserDefinedMaxwellExternalCurrent
100 !%Type block
101 !%Section Maxwell
102 !%Description
103 !%
104 !% Example:
105 !%
106 !% <tt>%UserDefinedMaxwellExternalCurrent
107 !% <br>&nbsp;&nbsp; current_parser | "expression_x_dir1" | "expression_y_dir1" | "expression_z_dir1"
108 !% <br>&nbsp;&nbsp; current_parser | "expression_x_dir2" | "expression_y_dir2" | "expression_z_dir2"
109 !% <br>&nbsp;&nbsp; current_td_function | "amplitude_j0_x" | "amplitude_j0_y" | "amplitude_j0_z" | omega |
110 !% envelope_td_function_name | phase
111 !% <br>%</tt>
112 !%
113 !% Description about UserDefinedMaxwellExternalCurrent follows
114 !%
115 !%Option current_parser 0
116 !% description follows
117 !%Option current_td_function 1
118 !% description follows
119 !%End
120
121 if (parse_block(namespace, 'UserDefinedMaxwellExternalCurrent', blk) == 0) then
122
123 ! find out how many lines (i.e. states) the block has
124 nlines = parse_block_n(blk)
125
126 st%external_current_number = nlines
127 safe_allocate(st%external_current_modus(1:nlines))
128 safe_allocate(st%external_current_string(1:space%dim, 1:nlines))
129 safe_allocate(st%external_current_amplitude(1:mesh%np, 1:space%dim, 1:nlines))
130 safe_allocate(st%external_current_td_function(1:nlines))
131 safe_allocate(st%external_current_omega(1:nlines))
132 safe_allocate(st%external_current_td_phase(1:nlines))
133
134 ! read all lines
135 do il = 1, nlines
136 ! Check that number of columns is four, five, six or seven.
137 ncols = parse_block_cols(blk, il - 1)
138 if ((ncols /= 4) .and. (ncols /= 5) .and. (ncols /= 6) .and. (ncols /= 7)) then
139 message(1) = 'Each line in the MaxwellExternalCurrent block must have'
140 message(2) = 'four, five, six or or seven columns.'
141 call messages_fatal(2, namespace=namespace)
142 end if
143
144 call parse_block_integer(blk, il - 1, 0, st%external_current_modus(il))
145
146 if (st%external_current_modus(il) == external_current_parser) then
147 ! parse formula string
148 do idir = 1, space%dim
149 call parse_block_string(blk, il - 1, idir, st%external_current_string(idir, il), convert_to_c=.true.)
150 end do
151 else if (st%external_current_modus(il) == external_current_td_function) then
152 do idir = 1, space%dim
153 call parse_block_string(blk, il - 1, idir, st%external_current_string(idir, il), convert_to_c=.true.)
154 do ip = 1, mesh%np
155 call mesh_r(mesh, ip, rr, coords = xx)
156 call parse_expression(j_vector(idir), dummy(idir), st%dim, xx, rr, m_zero, &
157 st%external_current_string(idir, il))
158 st%external_current_amplitude(ip, idir, il) = j_vector(idir)
159 end do
160 end do
161 call parse_block_float(blk, il-1, 4, omega, unit_one/units_inp%time)
162 st%external_current_omega(il) = omega
163 call parse_block_string(blk, il-1, 5, tdf_expression)
164 call tdf_read(st%external_current_td_function(il), namespace, trim(tdf_expression), ierr)
165 if (parse_block_cols(blk, il-1) > 6) then
166 call parse_block_string(blk, il-1, 6, phase_expression)
167 call tdf_read(st%external_current_td_phase(il), namespace, trim(phase_expression), ierr)
168 if (ierr /= 0) then
169 write(message(1),'(3A)') 'Error in the "', trim(tdf_expression), '" field defined in the TDExternalFields block:'
170 write(message(2),'(3A)') 'Time-dependent phase function "', trim(phase_expression), '" not found.'
171 call messages_warning(2, namespace=namespace)
172 end if
173 else
174 call tdf_init(st%external_current_td_phase(il))
175 end if
176 end if
177 end do
179 end if
180
181 call profiling_out('EXTERNAL_CURRENT_INIT')
182
183 pop_sub(external_current_init)
184 end subroutine external_current_init
185
186 !----------------------------------------------------------
187 subroutine external_current_calculation(st, space, mesh, time, current)
188 type(states_mxll_t), intent(in) :: st
189 class(space_t), intent(in) :: space
190 class(mesh_t), intent(in) :: mesh
191 real(real64), intent(in) :: time
192 real(real64), intent(inout) :: current(:,:)
193
194 integer :: ip, jn, idir
195 real(real64) :: xx(space%dim), rr, tt, j_vector(space%dim), dummy(space%dim), amp(space%dim)
196 complex(real64) :: exp_arg
197 real(real64) :: tmp_amp, phase
198
200
201 call profiling_in("EXTERNAL_CURRENT_CALC")
202
203 current(:,:) = m_zero
204 do jn = 1, st%external_current_number
205 if (st%external_current_modus(jn) == external_current_parser) then
206 do ip = 1, mesh%np
207 call mesh_r(mesh, ip, rr, coords = xx)
208 do idir = 1, space%dim
209 tt = time
210 call parse_expression(j_vector(idir), dummy(idir), space%dim, xx, rr, tt, &
211 trim(st%external_current_string(idir,jn)))
212 end do
213 current(ip, :) = current(ip, :) + j_vector
214 end do
215
216 else if (st%external_current_modus(jn) == external_current_td_function) then
217 exp_arg = st%external_current_omega(jn) * time + tdf(st%external_current_td_phase(jn),time)
218 phase = real(exp(-m_zi*exp_arg), real64)
219 tmp_amp = tdf(st%external_current_td_function(jn), time)
220 do ip = 1, mesh%np
221 amp = st%external_current_amplitude(ip, :, jn) * tmp_amp
222 current(ip, :) = current(ip, :) + amp(:) * phase
223 end do
224 end if
225 end do
226
227 call profiling_out("EXTERNAL_CURRENT_CALC")
228
230 end subroutine external_current_calculation
231
233
234!! Local Variables:
235!! mode: f90
236!! coding: utf-8
237!! End:
double exp(double __x) __attribute__((__nothrow__
subroutine, public get_rs_density_ext(st, space, mesh, time, rs_current_density_ext)
integer, parameter, public external_current_td_function
subroutine, public external_current_calculation(st, space, mesh, time, current)
subroutine, public external_current_init(st, space, namespace, mesh)
real(real64), parameter, public m_zero
Definition: global.F90:191
complex(real64), parameter, public m_zi
Definition: global.F90:205
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:341
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public parse_block_string(blk, l, c, res, convert_to_c)
Definition: parser.F90:810
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
subroutine, public build_rs_current_state(current_state, mesh, rs_current_state, ep_field, np)
subroutine, public tdf_init(f)
Definition: tdfunction.F90:390
subroutine, public tdf_read(f, namespace, function_name, ierr)
This function initializes "f" from the TDFunctions block.
Definition: tdfunction.F90:220
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
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
Definition: xc.F90:116
Describes mesh distribution to nodes.
Definition: mesh.F90:187
int true(void)