Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief subroutines for defining and creating Dirichlet type subdomains
10 : !> \par History
11 : !> 08.2014 created [Hossein Bani-Hashemian]
12 : !> 10.2015 completely revised [Hossein Bani-Hashemian]
13 : !> \author Mohammad Hossein Bani-Hashemian
14 : ! **************************************************************************************************
15 : MODULE dirichlet_bc_methods
16 :
17 : USE cp_log_handling, ONLY: cp_get_default_logger,&
18 : cp_logger_get_default_unit_nr,&
19 : cp_logger_type
20 : USE dirichlet_bc_types, ONLY: &
21 : AA_CUBOIDAL, AA_PLANAR, CIRCUMSCRIBED, CYLINDRICAL, INSCRIBED, PLANAR, &
22 : dbc_parameters_dealloc, dirichlet_bc_p_type, dirichlet_bc_type, tile_p_type, x_axis, &
23 : xy_plane, xz_plane, y_axis, yz_plane, z_axis
24 : USE kinds, ONLY: dp
25 : USE mathconstants, ONLY: pi,&
26 : twopi
27 : USE mathlib, ONLY: inv_3x3,&
28 : vector_product
29 : USE physcon, ONLY: angstrom
30 : USE ps_implicit_types, ONLY: MIXED_BC,&
31 : MIXED_PERIODIC_BC,&
32 : NEUMANN_BC,&
33 : PERIODIC_BC
34 : USE pw_grid_types, ONLY: pw_grid_type
35 : USE pw_methods, ONLY: pw_axpy,&
36 : pw_zero
37 : USE pw_poisson_types, ONLY: pw_poisson_parameter_type
38 : USE pw_pool_types, ONLY: pw_pool_type
39 : USE pw_types, ONLY: pw_r3d_rs_type
40 : USE rs_methods, ONLY: setup_grid_axes
41 : #include "../base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 : PRIVATE
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dirichlet_bc_methods'
46 :
47 : PUBLIC dirichlet_boundary_region_setup
48 :
49 : REAL(dp), PARAMETER, PRIVATE :: small_value = 1.0E-8_dp
50 : INTEGER, PARAMETER, PRIVATE :: FWROT = 1, &
51 : BWROT = -1
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Sets up the Dirichlet boundary condition
57 : !> \param pw_pool pool of plane wave grid
58 : !> \param poisson_params poisson_env parameters
59 : !> \param dbcs the DBC region to be created
60 : !> \par History
61 : !> 10.2014 created [Hossein Bani-Hashemian]
62 : !> \author Mohammad Hossein Bani-Hashemian
63 : ! **************************************************************************************************
64 54 : SUBROUTINE dirichlet_boundary_region_setup(pw_pool, poisson_params, dbcs)
65 :
66 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
67 : TYPE(pw_poisson_parameter_type), INTENT(INOUT) :: poisson_params
68 : TYPE(dirichlet_bc_p_type), ALLOCATABLE, &
69 : DIMENSION(:), INTENT(INOUT) :: dbcs
70 :
71 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dirichlet_boundary_region_setup'
72 :
73 : INTEGER :: apx_type, dbc_id, handle, ind_end, ind_start, j, l, n_aa_cuboidal, &
74 : n_aa_cylindrical, n_aa_planar, n_dbcs, n_planar, parallel_axis, parallel_plane, unit_nr
75 : INTEGER, DIMENSION(3) :: n_prtn, npts
76 54 : INTEGER, DIMENSION(:), POINTER :: aa_cylindrical_nsides
77 : LOGICAL :: is_periodic, verbose
78 : REAL(dp) :: base_radius, delta_alpha, frequency, &
79 : oscillating_fraction, phase, sigma, &
80 : thickness, v_D
81 54 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: x_glbl, x_locl, y_glbl, y_locl, z_glbl, &
82 54 : z_locl
83 : REAL(dp), DIMENSION(2) :: base_center
84 : REAL(dp), DIMENSION(2, 3) :: cell_xtnts
85 : REAL(dp), DIMENSION(3) :: dr
86 : TYPE(cp_logger_type), POINTER :: logger
87 : TYPE(pw_grid_type), POINTER :: pw_grid
88 :
89 54 : CALL timeset(routineN, handle)
90 :
91 54 : logger => cp_get_default_logger()
92 54 : IF (logger%para_env%is_source()) THEN
93 27 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
94 : ELSE
95 : unit_nr = -1
96 : END IF
97 :
98 216 : dr = pw_pool%pw_grid%dr
99 216 : npts = pw_pool%pw_grid%npts
100 54 : cell_xtnts = 0.0_dp
101 54 : cell_xtnts(2, 1) = npts(1)*dr(1)
102 54 : cell_xtnts(2, 2) = npts(2)*dr(2)
103 54 : cell_xtnts(2, 3) = npts(3)*dr(3)
104 :
105 54 : verbose = poisson_params%dbc_params%verbose_output
106 :
107 54 : n_aa_planar = poisson_params%dbc_params%n_aa_planar
108 54 : n_aa_cuboidal = poisson_params%dbc_params%n_aa_cuboidal
109 54 : n_planar = poisson_params%dbc_params%n_planar
110 54 : n_aa_cylindrical = poisson_params%dbc_params%n_aa_cylindrical
111 54 : aa_cylindrical_nsides => poisson_params%dbc_params%aa_cylindrical_nsides
112 54 : pw_grid => pw_pool%pw_grid
113 : CALL setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
114 66 : n_dbcs = n_aa_planar + n_aa_cuboidal + n_planar + SUM(aa_cylindrical_nsides)
115 :
116 92 : SELECT CASE (poisson_params%ps_implicit_params%boundary_condition)
117 : CASE (MIXED_BC, MIXED_PERIODIC_BC)
118 38 : IF (n_dbcs == 0) THEN
119 0 : CPABORT("No Dirichlet region is defined.")
120 : END IF
121 :
122 242 : ALLOCATE (dbcs(n_dbcs))
123 38 : IF (unit_nr > 0) THEN
124 19 : WRITE (unit_nr, '(/,T3,A,A,/,T3,A)') "POISSON| IMPLICIT (GENERALIZED) SOLVER ", REPEAT('-', 39), &
125 38 : "POISSON| Preparing Dirichlet regions ..."
126 : END IF
127 :
128 62 : DO j = 1, n_aa_planar
129 792 : ALLOCATE (dbcs(j)%dirichlet_bc)
130 96 : n_prtn = poisson_params%dbc_params%aa_planar_nprtn(:, j)
131 24 : dbc_id = AA_PLANAR + j
132 24 : v_D = poisson_params%dbc_params%aa_planar_vD(j)
133 24 : frequency = poisson_params%dbc_params%aa_planar_frequency(j)
134 24 : oscillating_fraction = poisson_params%dbc_params%aa_planar_osc_frac(j)
135 24 : phase = poisson_params%dbc_params%aa_planar_phase(j)
136 24 : sigma = poisson_params%dbc_params%aa_planar_sigma(j)
137 24 : thickness = poisson_params%dbc_params%aa_planar_thickness(j)
138 24 : parallel_plane = poisson_params%dbc_params%aa_planar_pplane(j)
139 24 : is_periodic = poisson_params%dbc_params%aa_planar_is_periodic(j)
140 :
141 24 : IF (unit_nr > 0) THEN
142 12 : WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", j
143 12 : WRITE (unit_nr, '(T3,A)') "POISSON| type : axis-aligned planar"
144 12 : WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_D, "[Eh/e]"
145 : END IF
146 : CALL aa_planar_dbc_setup(cell_xtnts, parallel_plane, &
147 : poisson_params%dbc_params%aa_planar_xxtnt(:, j), &
148 : poisson_params%dbc_params%aa_planar_yxtnt(:, j), &
149 : poisson_params%dbc_params%aa_planar_zxtnt(:, j), &
150 : thickness, sigma, v_D, oscillating_fraction, frequency, &
151 24 : phase, dbc_id, verbose, dbcs(j)%dirichlet_bc)
152 :
153 : CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
154 62 : x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
155 : END DO
156 :
157 38 : l = n_aa_planar
158 44 : DO j = l + 1, l + n_aa_cuboidal
159 198 : ALLOCATE (dbcs(j)%dirichlet_bc)
160 24 : n_prtn = poisson_params%dbc_params%aa_cuboidal_nprtn(:, j - l)
161 6 : dbc_id = AA_CUBOIDAL + j - l
162 6 : v_D = poisson_params%dbc_params%aa_cuboidal_vD(j - l)
163 6 : frequency = poisson_params%dbc_params%aa_cuboidal_frequency(j - l)
164 6 : oscillating_fraction = poisson_params%dbc_params%aa_cuboidal_osc_frac(j - l)
165 6 : phase = poisson_params%dbc_params%aa_cuboidal_phase(j - l)
166 6 : sigma = poisson_params%dbc_params%aa_cuboidal_sigma(j - l)
167 6 : is_periodic = poisson_params%dbc_params%aa_cuboidal_is_periodic(j - l)
168 :
169 6 : IF (unit_nr > 0) THEN
170 3 : WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", j
171 3 : WRITE (unit_nr, '(T3,A)') "POISSON| type : axis-aligned cuboidal"
172 3 : WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_D, "[Eh/e]"
173 : END IF
174 : CALL aa_cuboidal_dbc_setup(cell_xtnts, &
175 : poisson_params%dbc_params%aa_cuboidal_xxtnt(:, j - l), &
176 : poisson_params%dbc_params%aa_cuboidal_yxtnt(:, j - l), &
177 : poisson_params%dbc_params%aa_cuboidal_zxtnt(:, j - l), &
178 : sigma, v_D, oscillating_fraction, frequency, &
179 6 : phase, dbc_id, verbose, dbcs(j)%dirichlet_bc)
180 :
181 : CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
182 44 : x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
183 : END DO
184 :
185 38 : l = n_aa_planar + n_aa_cuboidal
186 44 : DO j = l + 1, l + n_planar
187 198 : ALLOCATE (dbcs(j)%dirichlet_bc)
188 24 : n_prtn = 1
189 18 : n_prtn(1:2) = poisson_params%dbc_params%planar_nprtn(:, j - l)
190 6 : dbc_id = PLANAR + j - l
191 6 : v_D = poisson_params%dbc_params%planar_vD(j - l)
192 6 : frequency = poisson_params%dbc_params%planar_frequency(j - l)
193 6 : oscillating_fraction = poisson_params%dbc_params%planar_osc_frac(j - l)
194 6 : phase = poisson_params%dbc_params%planar_phase(j - l)
195 6 : sigma = poisson_params%dbc_params%planar_sigma(j - l)
196 6 : thickness = poisson_params%dbc_params%planar_thickness(j - l)
197 6 : is_periodic = poisson_params%dbc_params%planar_is_periodic(j - l)
198 :
199 6 : IF (unit_nr > 0) THEN
200 3 : WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", j
201 3 : WRITE (unit_nr, '(T3,A)') "POISSON| type : planar"
202 3 : WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_D, "[Eh/e]"
203 : END IF
204 : CALL arbitrary_planar_dbc_setup(cell_xtnts, &
205 : poisson_params%dbc_params%planar_Avtx(:, j - l), &
206 : poisson_params%dbc_params%planar_Bvtx(:, j - l), &
207 : poisson_params%dbc_params%planar_Cvtx(:, j - l), &
208 : thickness, sigma, v_D, oscillating_fraction, frequency, &
209 6 : phase, dbc_id, verbose, dbcs(j)%dirichlet_bc)
210 :
211 : CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
212 44 : x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
213 : END DO
214 :
215 38 : l = n_aa_planar + n_aa_cuboidal + n_planar
216 104 : DO j = 1, n_aa_cylindrical
217 12 : ind_start = l + 1
218 12 : ind_end = l + aa_cylindrical_nsides(j)
219 :
220 48 : n_prtn = 1
221 36 : n_prtn(1:2) = poisson_params%dbc_params%aa_cylindrical_nprtn(:, j)
222 12 : parallel_axis = poisson_params%dbc_params%aa_cylindrical_paxis(j)
223 12 : apx_type = poisson_params%dbc_params%aa_cylindrical_apxtyp(j)
224 36 : base_center = poisson_params%dbc_params%aa_cylindrical_bctr(:, j)
225 12 : base_radius = poisson_params%dbc_params%aa_cylindrical_brad(j)
226 12 : thickness = poisson_params%dbc_params%aa_cylindrical_thickness(j)
227 12 : delta_alpha = poisson_params%dbc_params%aa_cylindrical_sgap(j)
228 12 : sigma = poisson_params%dbc_params%aa_cylindrical_sigma(j)
229 12 : v_D = poisson_params%dbc_params%aa_cylindrical_vD(j)
230 12 : frequency = poisson_params%dbc_params%aa_cylindrical_frequency(j)
231 12 : oscillating_fraction = poisson_params%dbc_params%aa_cylindrical_osc_frac(j)
232 12 : phase = poisson_params%dbc_params%aa_cylindrical_phase(j)
233 12 : is_periodic = poisson_params%dbc_params%aa_cylindrical_is_periodic(j)
234 :
235 12 : IF (unit_nr > 0) THEN
236 6 : WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", l + j
237 6 : WRITE (unit_nr, '(T3,A)') "POISSON| type : axis-aligned cylindrical"
238 6 : WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_D, "[Eh/e]"
239 : END IF
240 : CALL aa_cylindrical_dbc_setup(pw_pool, cell_xtnts, x_locl, y_locl, z_locl, apx_type, parallel_axis, &
241 : poisson_params%dbc_params%aa_cylindrical_xtnt(:, j), &
242 : base_center, base_radius, thickness, delta_alpha, sigma, v_D, &
243 : oscillating_fraction, frequency, phase, &
244 12 : n_prtn, is_periodic, verbose, dbcs(ind_start:ind_end))
245 :
246 50 : l = l + aa_cylindrical_nsides(j)
247 : END DO
248 :
249 : CASE (PERIODIC_BC, NEUMANN_BC)
250 : ! do nothing
251 : END SELECT
252 :
253 : ! we won't need parameters anymore so deallocate them
254 54 : CALL dbc_parameters_dealloc(poisson_params%dbc_params)
255 :
256 54 : CALL timestop(handle)
257 :
258 108 : END SUBROUTINE dirichlet_boundary_region_setup
259 :
260 : ! **************************************************************************************************
261 : !> \brief Partitions a dirichlet_bc_type into tiles
262 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
263 : !> \param n_prtn vector of size 3 specifying the number of times that the x, y and
264 : !> z interval (defining the region) should be partitioned into
265 : !> \param pw_pool pool of plane-wave grid
266 : !> \param cell_xtnts extents of the simulation cell
267 : !> \param x_locl x grid vector of the simulation box local to this process
268 : !> \param y_locl y grid vector of the simulation box local to this process
269 : !> \param z_locl z grid vector of the simulation box local to this process
270 : !> \param is_periodic whether or not the tile is periodic with a periodicity of one unit cell
271 : !> \param verbose whether or not to print out the coordinates of the vertices
272 : !> \param dirichlet_bc the dirichlet_bc object to be partitioned
273 : !> \par History
274 : !> 10.2014 created [Hossein Bani-Hashemian]
275 : !> \author Mohammad Hossein Bani-Hashemian
276 : ! **************************************************************************************************
277 128 : SUBROUTINE dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
278 : x_locl, y_locl, z_locl, is_periodic, verbose, dirichlet_bc)
279 :
280 : REAL(dp), INTENT(IN) :: sigma
281 : INTEGER, DIMENSION(3), INTENT(IN) :: n_prtn
282 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
283 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
284 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
285 : LOGICAL, INTENT(IN) :: is_periodic, verbose
286 : TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
287 :
288 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dirichlet_bc_partition'
289 :
290 : INTEGER :: handle, i, k, n_tiles, unit_nr
291 : REAL(dp) :: cyl_maxval, tile_volume
292 128 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: tile_maxvals
293 : TYPE(cp_logger_type), POINTER :: logger
294 : TYPE(pw_r3d_rs_type), POINTER :: cylinder_pw
295 128 : TYPE(tile_p_type), DIMENSION(:), POINTER :: tiles
296 :
297 128 : CALL timeset(routineN, handle)
298 :
299 128 : logger => cp_get_default_logger()
300 128 : IF (logger%para_env%is_source()) THEN
301 64 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
302 : ELSE
303 : unit_nr = -1
304 : END IF
305 :
306 134 : SELECT CASE (dirichlet_bc%dbc_geom)
307 : CASE (PLANAR)
308 :
309 6 : CALL arbitrary_dbc_partition(dirichlet_bc%vertices, n_prtn(1:2), tiles)
310 6 : n_tiles = SIZE(tiles)
311 6 : dirichlet_bc%n_tiles = n_tiles
312 36 : ALLOCATE (dirichlet_bc%tiles(n_tiles))
313 42 : dirichlet_bc%tiles = tiles
314 6 : IF ((unit_nr > 0) .AND. verbose) WRITE (unit_nr, '(T7,A,I5)') "number of partitions : ", n_tiles
315 :
316 12 : ALLOCATE (tile_maxvals(n_tiles))
317 24 : tile_maxvals = 0.0_dp
318 24 : DO k = 1, n_tiles
319 18 : ALLOCATE (dirichlet_bc%tiles(k)%tile%tile_pw)
320 18 : CALL pw_pool%create_pw(dirichlet_bc%tiles(k)%tile%tile_pw)
321 18 : CALL pw_zero(dirichlet_bc%tiles(k)%tile%tile_pw)
322 :
323 18 : IF ((unit_nr > 0) .AND. verbose) THEN
324 0 : WRITE (unit_nr, '(T7,A,I5)') "tile", k
325 0 : DO i = 1, 8
326 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') &
327 0 : " vertex ", i, angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
328 : END DO
329 : END IF
330 :
331 : CALL arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, dirichlet_bc%tiles(k)%tile%vertices, &
332 18 : sigma, dirichlet_bc%tiles(k)%tile%tile_pw)
333 2270862 : tile_volume = SUM(dirichlet_bc%tiles(k)%tile%tile_pw%array)
334 18 : CALL pw_pool%pw_grid%para%group%sum(tile_volume)
335 24 : dirichlet_bc%tiles(k)%tile%volume = tile_volume
336 : END DO
337 :
338 6 : IF ((unit_nr > 0) .AND. verbose) WRITE (unit_nr, '(T3,A)') REPEAT('=', 78)
339 :
340 6 : DEALLOCATE (tiles)
341 :
342 : CASE (CYLINDRICAL)
343 :
344 92 : CALL arbitrary_dbc_partition(dirichlet_bc%vertices, n_prtn(1:2), tiles)
345 92 : n_tiles = SIZE(tiles)
346 92 : dirichlet_bc%n_tiles = n_tiles
347 418 : ALLOCATE (dirichlet_bc%tiles(n_tiles))
348 376 : dirichlet_bc%tiles = tiles
349 92 : IF ((unit_nr > 0) .AND. verbose) WRITE (unit_nr, '(T7,A,I5)') "number of partitions : ", n_tiles
350 :
351 92 : NULLIFY (cylinder_pw)
352 92 : ALLOCATE (cylinder_pw)
353 92 : CALL pw_pool%create_pw(cylinder_pw)
354 92 : CALL pw_zero(cylinder_pw)
355 :
356 234 : DO k = 1, n_tiles
357 142 : ALLOCATE (dirichlet_bc%tiles(k)%tile%tile_pw)
358 142 : CALL pw_pool%create_pw(dirichlet_bc%tiles(k)%tile%tile_pw)
359 142 : CALL pw_zero(dirichlet_bc%tiles(k)%tile%tile_pw)
360 :
361 142 : IF ((unit_nr > 0) .AND. verbose) THEN
362 0 : WRITE (unit_nr, '(T7,A,I5)') "tile", k
363 0 : DO i = 1, 8
364 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') &
365 0 : " vertex ", i, angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
366 : END DO
367 : END IF
368 :
369 : CALL arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, dirichlet_bc%tiles(k)%tile%vertices, &
370 142 : sigma, dirichlet_bc%tiles(k)%tile%tile_pw)
371 142 : CALL pw_axpy(dirichlet_bc%tiles(k)%tile%tile_pw, cylinder_pw)
372 13544386 : tile_volume = SUM(dirichlet_bc%tiles(k)%tile%tile_pw%array)
373 142 : CALL pw_pool%pw_grid%para%group%sum(tile_volume)
374 234 : dirichlet_bc%tiles(k)%tile%volume = tile_volume
375 : END DO
376 :
377 7255676 : cyl_maxval = MAXVAL(cylinder_pw%array)
378 92 : CALL pw_pool%pw_grid%para%group%max(cyl_maxval)
379 92 : IF (cyl_maxval > 1.01_dp) THEN
380 : CALL cp_abort(__LOCATION__, &
381 : "Inaccurate approximation of unit step function at cylindrical Dirichlet region. "// &
382 : "Increase the number of sides (N_SIDES) and/or the base radius. Alternatively, "// &
383 0 : "one can increase DELTA_ALPHA (see the reference manual).")
384 : END IF
385 :
386 92 : CALL pw_pool%give_back_pw(cylinder_pw)
387 92 : DEALLOCATE (cylinder_pw)
388 :
389 92 : IF ((unit_nr > 0) .AND. verbose) WRITE (unit_nr, '(T3,A)') REPEAT('=', 78)
390 :
391 92 : DEALLOCATE (tiles)
392 :
393 : CASE (AA_PLANAR, AA_CUBOIDAL)
394 :
395 30 : CALL aa_dbc_partition(dirichlet_bc%vertices, n_prtn, tiles)
396 30 : n_tiles = SIZE(tiles)
397 30 : dirichlet_bc%n_tiles = n_tiles
398 154 : ALLOCATE (dirichlet_bc%tiles(n_tiles))
399 158 : dirichlet_bc%tiles = tiles
400 30 : IF ((unit_nr > 0) .AND. verbose) WRITE (unit_nr, '(T7,A,I5)') "number of partitions : ", n_tiles
401 :
402 94 : DO k = 1, n_tiles
403 64 : ALLOCATE (dirichlet_bc%tiles(k)%tile%tile_pw)
404 64 : CALL pw_pool%create_pw(dirichlet_bc%tiles(k)%tile%tile_pw)
405 64 : CALL pw_zero(dirichlet_bc%tiles(k)%tile%tile_pw)
406 :
407 64 : IF ((unit_nr > 0) .AND. verbose) THEN
408 0 : WRITE (unit_nr, '(T7,A,I5)') "tile", k
409 0 : DO i = 1, 8
410 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') &
411 0 : " vertex ", i, angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
412 : END DO
413 : END IF
414 :
415 : CALL aa_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, dirichlet_bc%tiles(k)%tile%vertices, &
416 64 : sigma, is_periodic, dirichlet_bc%tiles(k)%tile%tile_pw)
417 8269897 : tile_volume = SUM(dirichlet_bc%tiles(k)%tile%tile_pw%array)
418 64 : CALL pw_pool%pw_grid%para%group%sum(tile_volume)
419 94 : dirichlet_bc%tiles(k)%tile%volume = tile_volume
420 : END DO
421 :
422 30 : IF ((unit_nr > 0) .AND. verbose) WRITE (unit_nr, '(T3,A)') REPEAT('=', 78)
423 :
424 158 : DEALLOCATE (tiles)
425 :
426 : END SELECT
427 :
428 128 : CALL timestop(handle)
429 :
430 256 : END SUBROUTINE dirichlet_bc_partition
431 :
432 : ! **************************************************************************************************
433 : !> \brief Creates an axis-aligned planar Dirichlet region.
434 : !> \param cell_xtnts extents of the simulation cell
435 : !> \param parallel_plane the coordinate plane that the region is parallel to
436 : !> \param x_xtnt the x extent of the planar region
437 : !> \param y_xtnt the y extent of the planar region
438 : !> \param z_xtnt the z extent of the planar region
439 : !> \param thickness the thickness of the region
440 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
441 : !> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
442 : !> \param osc_frac ...
443 : !> \param frequency ...
444 : !> \param phase ...
445 : !> \param dbc_id unique ID for the Dirichlet region
446 : !> \param verbose whether or not to print out the coordinates of the vertices
447 : !> \param dirichlet_bc the dirichlet_bc object to be created
448 : !> \par History
449 : !> 08.2014 created [Hossein Bani-Hashemian]
450 : !> \author Mohammad Hossein Bani-Hashemian
451 : ! **************************************************************************************************
452 24 : SUBROUTINE aa_planar_dbc_setup(cell_xtnts, parallel_plane, x_xtnt, y_xtnt, z_xtnt, &
453 : thickness, sigma, v_D, osc_frac, frequency, phase, dbc_id, verbose, dirichlet_bc)
454 :
455 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
456 : INTEGER, INTENT(IN) :: parallel_plane
457 : REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, y_xtnt, z_xtnt
458 : REAL(dp), INTENT(IN) :: thickness, sigma, v_D, osc_frac, &
459 : frequency, phase
460 : INTEGER, INTENT(IN) :: dbc_id
461 : LOGICAL, INTENT(IN) :: verbose
462 : TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
463 :
464 : CHARACTER(LEN=*), PARAMETER :: routineN = 'aa_planar_dbc_setup'
465 : LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .TRUE.
466 :
467 : INTEGER :: handle, i, n_forb_xtnts, unit_nr
468 : LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
469 : forb_xtnt4, forb_xtnt5, forb_xtnt6
470 : REAL(dp) :: xlb, xub, ylb, yub, zlb, zub
471 : TYPE(cp_logger_type), POINTER :: logger
472 :
473 24 : CALL timeset(routineN, handle)
474 :
475 24 : logger => cp_get_default_logger()
476 24 : IF (logger%para_env%is_source()) THEN
477 12 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
478 : ELSE
479 : unit_nr = -1
480 : END IF
481 :
482 24 : xlb = x_xtnt(1); xub = x_xtnt(2)
483 24 : ylb = y_xtnt(1); yub = y_xtnt(2)
484 24 : zlb = z_xtnt(1); zub = z_xtnt(2)
485 :
486 26 : SELECT CASE (parallel_plane)
487 : CASE (xy_plane)
488 2 : zlb = zlb - thickness*0.5_dp
489 2 : zub = zub + thickness*0.5_dp
490 : CASE (xz_plane)
491 2 : ylb = ylb - thickness*0.5_dp
492 2 : yub = yub + thickness*0.5_dp
493 : CASE (yz_plane)
494 20 : xlb = xlb - thickness*0.5_dp
495 24 : xub = xub + thickness*0.5_dp
496 : END SELECT
497 :
498 24 : forb_xtnt1 = xlb < cell_xtnts(1, 1)
499 24 : forb_xtnt2 = xub > cell_xtnts(2, 1)
500 24 : forb_xtnt3 = ylb < cell_xtnts(1, 2)
501 24 : forb_xtnt4 = yub > cell_xtnts(2, 2)
502 24 : forb_xtnt5 = zlb < cell_xtnts(1, 3)
503 24 : forb_xtnt6 = zub > cell_xtnts(2, 3)
504 : n_forb_xtnts = COUNT([forb_xtnt1, forb_xtnt2, forb_xtnt3, &
505 168 : forb_xtnt4, forb_xtnt5, forb_xtnt6] .EQV. test_forb_xtnts)
506 24 : IF (n_forb_xtnts > 0) THEN
507 : CALL cp_abort(__LOCATION__, &
508 : "The given extents for the Dirichlet region are outside the range of "// &
509 0 : "the simulation cell.")
510 : END IF
511 :
512 24 : dirichlet_bc%v_D = v_D
513 24 : dirichlet_bc%osc_frac = osc_frac
514 24 : dirichlet_bc%frequency = frequency
515 24 : dirichlet_bc%phase = phase
516 24 : dirichlet_bc%dbc_id = dbc_id
517 24 : dirichlet_bc%dbc_geom = AA_PLANAR
518 24 : dirichlet_bc%smoothing_width = sigma
519 24 : dirichlet_bc%n_tiles = 1
520 :
521 96 : dirichlet_bc%vertices(1:3, 1) = [xlb, ylb, zlb]
522 96 : dirichlet_bc%vertices(1:3, 2) = [xub, ylb, zlb]
523 96 : dirichlet_bc%vertices(1:3, 3) = [xub, yub, zlb]
524 96 : dirichlet_bc%vertices(1:3, 4) = [xlb, yub, zlb]
525 96 : dirichlet_bc%vertices(1:3, 5) = [xlb, yub, zub]
526 96 : dirichlet_bc%vertices(1:3, 6) = [xlb, ylb, zub]
527 96 : dirichlet_bc%vertices(1:3, 7) = [xub, ylb, zub]
528 96 : dirichlet_bc%vertices(1:3, 8) = [xub, yub, zub]
529 :
530 24 : IF ((unit_nr > 0) .AND. verbose) THEN
531 0 : WRITE (unit_nr, '(T3,A,A)') "======== verbose ", REPEAT('=', 61)
532 0 : DO i = 1, 8
533 0 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, angstrom*dirichlet_bc%vertices(:, i)
534 : END DO
535 : END IF
536 :
537 24 : CALL timestop(handle)
538 :
539 24 : END SUBROUTINE aa_planar_dbc_setup
540 :
541 : ! **************************************************************************************************
542 : !> \brief Creates an axis-aligned cuboidal Dirichlet region.
543 : !> \param cell_xtnts extents of the simulation cell
544 : !> \param x_xtnt the x extent of the cuboidal region
545 : !> \param y_xtnt the y extent of the cuboidal region
546 : !> \param z_xtnt the z extent of the cuboidal region
547 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
548 : !> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
549 : !> \param osc_frac ...
550 : !> \param frequency ...
551 : !> \param phase ...
552 : !> \param dbc_id unique ID for the Dirichlet region
553 : !> \param verbose whether or not to print out the coordinates of the vertices
554 : !> \param dirichlet_bc the dirichlet_bc object to be created
555 : !> \par History
556 : !> 12.2014 created [Hossein Bani-Hashemian]
557 : !> \author Mohammad Hossein Bani-Hashemian
558 : ! **************************************************************************************************
559 6 : SUBROUTINE aa_cuboidal_dbc_setup(cell_xtnts, x_xtnt, y_xtnt, z_xtnt, &
560 : sigma, v_D, osc_frac, frequency, phase, dbc_id, verbose, dirichlet_bc)
561 :
562 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
563 : REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, y_xtnt, z_xtnt
564 : REAL(dp), INTENT(IN) :: sigma, v_D, osc_frac, frequency, phase
565 : INTEGER, INTENT(IN) :: dbc_id
566 : LOGICAL, INTENT(IN) :: verbose
567 : TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
568 :
569 : CHARACTER(LEN=*), PARAMETER :: routineN = 'aa_cuboidal_dbc_setup'
570 : LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .TRUE.
571 :
572 : INTEGER :: handle, i, n_forb_xtnts, unit_nr
573 : LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
574 : forb_xtnt4, forb_xtnt5, forb_xtnt6
575 : TYPE(cp_logger_type), POINTER :: logger
576 :
577 6 : CALL timeset(routineN, handle)
578 :
579 6 : logger => cp_get_default_logger()
580 6 : IF (logger%para_env%is_source()) THEN
581 3 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
582 : ELSE
583 : unit_nr = -1
584 : END IF
585 :
586 6 : forb_xtnt1 = x_xtnt(1) < cell_xtnts(1, 1)
587 6 : forb_xtnt2 = x_xtnt(2) > cell_xtnts(2, 1)
588 6 : forb_xtnt3 = y_xtnt(1) < cell_xtnts(1, 2)
589 6 : forb_xtnt4 = y_xtnt(2) > cell_xtnts(2, 2)
590 6 : forb_xtnt5 = z_xtnt(1) < cell_xtnts(1, 3)
591 6 : forb_xtnt6 = z_xtnt(2) > cell_xtnts(2, 3)
592 : n_forb_xtnts = COUNT([forb_xtnt1, forb_xtnt2, forb_xtnt3, &
593 42 : forb_xtnt4, forb_xtnt5, forb_xtnt6] .EQV. test_forb_xtnts)
594 6 : IF (n_forb_xtnts > 0) THEN
595 : CALL cp_abort(__LOCATION__, &
596 : "The given extents for the Dirichlet region are outside the range of "// &
597 0 : "the simulation cell.")
598 : END IF
599 :
600 6 : dirichlet_bc%v_D = v_D
601 6 : dirichlet_bc%osc_frac = osc_frac
602 6 : dirichlet_bc%frequency = frequency
603 6 : dirichlet_bc%phase = phase
604 6 : dirichlet_bc%dbc_id = dbc_id
605 6 : dirichlet_bc%dbc_geom = AA_CUBOIDAL
606 6 : dirichlet_bc%smoothing_width = sigma
607 6 : dirichlet_bc%n_tiles = 1
608 :
609 24 : dirichlet_bc%vertices(1:3, 1) = [x_xtnt(1), y_xtnt(1), z_xtnt(1)]
610 24 : dirichlet_bc%vertices(1:3, 2) = [x_xtnt(2), y_xtnt(1), z_xtnt(1)]
611 24 : dirichlet_bc%vertices(1:3, 3) = [x_xtnt(2), y_xtnt(2), z_xtnt(1)]
612 24 : dirichlet_bc%vertices(1:3, 4) = [x_xtnt(1), y_xtnt(2), z_xtnt(1)]
613 24 : dirichlet_bc%vertices(1:3, 5) = [x_xtnt(1), y_xtnt(2), z_xtnt(2)]
614 24 : dirichlet_bc%vertices(1:3, 6) = [x_xtnt(1), y_xtnt(1), z_xtnt(2)]
615 24 : dirichlet_bc%vertices(1:3, 7) = [x_xtnt(2), y_xtnt(1), z_xtnt(2)]
616 24 : dirichlet_bc%vertices(1:3, 8) = [x_xtnt(2), y_xtnt(2), z_xtnt(2)]
617 :
618 6 : IF ((unit_nr > 0) .AND. verbose) THEN
619 0 : WRITE (unit_nr, '(T3,A,A)') "======== verbose ", REPEAT('=', 61)
620 0 : DO i = 1, 8
621 0 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, angstrom*dirichlet_bc%vertices(:, i)
622 : END DO
623 : END IF
624 :
625 6 : CALL timestop(handle)
626 :
627 6 : END SUBROUTINE aa_cuboidal_dbc_setup
628 :
629 : ! **************************************************************************************************
630 : !> \brief Creates an arbitrary (rectangular-shaped) planar Dirichlet region given
631 : !> the coordinates of its vertices.
632 : !> \param cell_xtnts extents of the simulation cell
633 : !> \param A coordinates of the vertex A
634 : !> \param B coordinates of the vertex B
635 : !> \param C coordinates of the vertex C
636 : !> \param thickness the thickness of the region
637 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
638 : !> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
639 : !> \param osc_frac ...
640 : !> \param frequency ...
641 : !> \param phase ...
642 : !> \param dbc_id unique ID for the Dirichlet region
643 : !> \param verbose whether or not to print out the coordinates of the vertices
644 : !> \param dirichlet_bc the dirichlet_bc object to be created
645 : !> \par History
646 : !> 08.2014 created [Hossein Bani-Hashemian]
647 : !> \author Mohammad Hossein Bani-Hashemian
648 : ! **************************************************************************************************
649 6 : SUBROUTINE arbitrary_planar_dbc_setup(cell_xtnts, A, B, C, thickness, &
650 : sigma, v_D, osc_frac, frequency, phase, dbc_id, verbose, dirichlet_bc)
651 :
652 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
653 : REAL(dp), DIMENSION(3), INTENT(IN) :: A, B, C
654 : REAL(dp), INTENT(IN) :: thickness, sigma, v_D, osc_frac, &
655 : frequency, phase
656 : INTEGER, INTENT(IN) :: dbc_id
657 : LOGICAL, INTENT(IN) :: verbose
658 : TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
659 :
660 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arbitrary_planar_dbc_setup'
661 :
662 : INTEGER :: handle, i, unit_nr
663 : LOGICAL :: A_is_inside, are_coplanar, &
664 : are_orthogonal, B_is_inside, &
665 : C_is_inside, D_is_inside, is_rectangle
666 : REAL(dp) :: dist1, dist2, dist3, dist4
667 : REAL(dp), DIMENSION(3) :: AB, AC, AD, BC, cm, D, normal_vector, &
668 : unit_normal
669 : TYPE(cp_logger_type), POINTER :: logger
670 :
671 6 : CALL timeset(routineN, handle)
672 :
673 6 : logger => cp_get_default_logger()
674 6 : IF (logger%para_env%is_source()) THEN
675 3 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
676 : ELSE
677 : unit_nr = -1
678 : END IF
679 :
680 24 : D = A + (C - B)
681 78 : AB = B - A; AC = C - A; AD = D - A; BC = C - B
682 60 : are_orthogonal = ABS(DOT_PRODUCT(AB, BC)/(SQRT(SUM(AB**2))*SQRT(SUM(BC**2)))) <= small_value
683 6 : IF (.NOT. are_orthogonal) THEN
684 : CALL cp_abort(__LOCATION__, "The given vertices for defining a "// &
685 0 : "planar Dirichlet region do not form orthogonal edges.")
686 : END IF
687 :
688 6 : normal_vector = vector_product(AB, AC)
689 42 : unit_normal = normal_vector/SQRT(SUM(normal_vector**2))
690 :
691 : A_is_inside = (A(1) > cell_xtnts(1, 1) .AND. A(1) < cell_xtnts(2, 1)) .AND. &
692 : (A(2) > cell_xtnts(1, 2) .AND. A(2) < cell_xtnts(2, 2)) .AND. &
693 6 : (A(3) > cell_xtnts(1, 3) .AND. A(3) < cell_xtnts(2, 3))
694 : B_is_inside = (B(1) > cell_xtnts(1, 1) .AND. B(1) < cell_xtnts(2, 1)) .AND. &
695 : (B(2) > cell_xtnts(1, 2) .AND. B(2) < cell_xtnts(2, 2)) .AND. &
696 6 : (B(3) > cell_xtnts(1, 3) .AND. B(3) < cell_xtnts(2, 3))
697 : C_is_inside = (C(1) > cell_xtnts(1, 1) .AND. C(1) < cell_xtnts(2, 1)) .AND. &
698 : (C(2) > cell_xtnts(1, 2) .AND. C(2) < cell_xtnts(2, 2)) .AND. &
699 6 : (C(3) > cell_xtnts(1, 3) .AND. C(3) < cell_xtnts(2, 3))
700 : D_is_inside = (D(1) > cell_xtnts(1, 1) .AND. D(1) < cell_xtnts(2, 1)) .AND. &
701 : (D(2) > cell_xtnts(1, 2) .AND. D(2) < cell_xtnts(2, 2)) .AND. &
702 6 : (D(3) > cell_xtnts(1, 3) .AND. D(3) < cell_xtnts(2, 3))
703 6 : IF (.NOT. (A_is_inside .AND. B_is_inside .AND. C_is_inside .AND. D_is_inside)) THEN
704 : CALL cp_abort(__LOCATION__, &
705 : "At least one of the given vertices for defining a planar Dirichlet "// &
706 0 : "region is outside the simulation box.")
707 : END IF
708 :
709 24 : are_coplanar = ABS(DOT_PRODUCT(vector_product(AB, AC), AD)) <= small_value
710 6 : IF (.NOT. are_coplanar) THEN
711 0 : CPABORT("The given vertices for defining a planar Dirichlet region are not coplanar.")
712 : END IF
713 :
714 24 : cm = (A + B + C + D)/4.0_dp
715 24 : dist1 = SUM((A - cm)**2)
716 24 : dist2 = SUM((B - cm)**2)
717 24 : dist3 = SUM((C - cm)**2)
718 24 : dist4 = SUM((D - cm)**2)
719 : is_rectangle = (ABS(dist1 - dist2) <= small_value) .AND. &
720 : (ABS(dist1 - dist3) <= small_value) .AND. &
721 6 : (ABS(dist1 - dist4) <= small_value)
722 : IF (.NOT. is_rectangle) THEN
723 : CALL cp_abort(__LOCATION__, "The given vertices for defining "// &
724 0 : "a planar Dirichlet region do not form a rectangle.")
725 : END IF
726 :
727 6 : dirichlet_bc%v_D = v_D
728 6 : dirichlet_bc%osc_frac = osc_frac
729 6 : dirichlet_bc%frequency = frequency
730 6 : dirichlet_bc%phase = phase
731 6 : dirichlet_bc%dbc_id = dbc_id
732 6 : dirichlet_bc%dbc_geom = PLANAR
733 6 : dirichlet_bc%smoothing_width = sigma
734 6 : dirichlet_bc%n_tiles = 1
735 :
736 24 : dirichlet_bc%vertices(1:3, 1) = A - 0.5_dp*thickness*unit_normal
737 24 : dirichlet_bc%vertices(1:3, 2) = B - 0.5_dp*thickness*unit_normal
738 24 : dirichlet_bc%vertices(1:3, 3) = C - 0.5_dp*thickness*unit_normal
739 24 : dirichlet_bc%vertices(1:3, 4) = D - 0.5_dp*thickness*unit_normal
740 24 : dirichlet_bc%vertices(1:3, 5) = D + 0.5_dp*thickness*unit_normal
741 24 : dirichlet_bc%vertices(1:3, 6) = A + 0.5_dp*thickness*unit_normal
742 24 : dirichlet_bc%vertices(1:3, 7) = B + 0.5_dp*thickness*unit_normal
743 24 : dirichlet_bc%vertices(1:3, 8) = C + 0.5_dp*thickness*unit_normal
744 :
745 6 : IF ((unit_nr > 0) .AND. verbose) THEN
746 0 : WRITE (unit_nr, '(T3,A,A)') "======== verbose ", REPEAT('=', 61)
747 0 : DO i = 1, 8
748 0 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, angstrom*dirichlet_bc%vertices(:, i)
749 : END DO
750 : END IF
751 :
752 6 : CALL timestop(handle)
753 :
754 6 : END SUBROUTINE arbitrary_planar_dbc_setup
755 :
756 : ! **************************************************************************************************
757 : !> \brief Creates an x-axis-aligned cylindrical Dirichlet region. The cylindrical
758 : !> shape is approximated by an n-gonal (circumscribed or inscribed) uniform prism-
759 : !> like object. The circular base is approximated by a polygon. A second polygon is
760 : !> then created by rotating the first one by delta_alpha radians. Let's name the
761 : !> vertices of the first and the second polygon as v1, v2, v3, ... vn-1, vn and
762 : !> w1, w2, w3, ... wn-1, wn, respectively. Then w1v2, w2v3, ... wn-1vn, wnv1 form
763 : !> the edges of the cross-section of the approximating object. This way the dbcs
764 : !> do not share any edge.
765 : !> \param pw_pool pool of plane wave grid
766 : !> \param cell_xtnts extents of the simulation cell
767 : !> \param x_locl x grid vector of the simulation box local to this process
768 : !> \param y_locl y grid vector of the simulation box local to this process
769 : !> \param z_locl z grid vector of the simulation box local to this process
770 : !> \param apx_type the type of the n-gonal prism approximating the cylinder
771 : !> \param parallel_axis the coordinate axis that the cylindrical region extends along
772 : !> \param xtnt the x extent of the cylindrical region along its central axis
773 : !> \param base_center the y and z coordinates of the cylinder's base
774 : !> \param base_radius the radius of the cylinder's base
775 : !> \param thickness the thickness of the region
776 : !> \param delta_alpha ...
777 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
778 : !> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
779 : !> \param osc_frac ...
780 : !> \param frequency ...
781 : !> \param phase ...
782 : !> \param n_prtn number of partitions along the edges of the faces of the prism
783 : !> \param is_periodic whether or not the tile is periodic with a periodicity of one unit cell
784 : !> \param verbose whether or not to print out the coordinates of the vertices
785 : !> \param dbcs the x-axis-aligned cylindrical gate region to be created
786 : !> \par History
787 : !> 08.2014 created [Hossein Bani-Hashemian]
788 : !> \author Mohammad Hossein Bani-Hashemian
789 : ! **************************************************************************************************
790 12 : SUBROUTINE aa_cylindrical_dbc_setup(pw_pool, cell_xtnts, x_locl, y_locl, z_locl, apx_type, parallel_axis, &
791 : xtnt, base_center, base_radius, thickness, delta_alpha, &
792 12 : sigma, v_D, osc_frac, frequency, phase, n_prtn, is_periodic, verbose, dbcs)
793 :
794 : TYPE(pw_pool_type), POINTER :: pw_pool
795 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
796 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
797 : INTEGER, INTENT(IN), OPTIONAL :: apx_type
798 : INTEGER, INTENT(IN) :: parallel_axis
799 : REAL(dp), DIMENSION(2), INTENT(IN) :: xtnt, base_center
800 : REAL(dp), INTENT(IN) :: base_radius, thickness, delta_alpha, &
801 : sigma, v_D, osc_frac, frequency, phase
802 : INTEGER, DIMENSION(3), INTENT(IN) :: n_prtn
803 : LOGICAL, INTENT(IN) :: is_periodic, verbose
804 : TYPE(dirichlet_bc_p_type), DIMENSION(:), &
805 : INTENT(INOUT) :: dbcs
806 :
807 : CHARACTER(LEN=*), PARAMETER :: routineN = 'aa_cylindrical_dbc_setup'
808 :
809 : INTEGER :: handle, i, intern_apx_type, j, n_dbcs, &
810 : unit_nr
811 : LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
812 : forb_xtnt4
813 : REAL(dp) :: cntrlaxis_lb, cntrlaxis_ub, h, Lx, Ly, &
814 : Lz, theta
815 : REAL(dp), DIMENSION(3) :: A, B, C, D, normal_vec, unit_normal
816 : TYPE(cp_logger_type), POINTER :: logger
817 24 : REAL(dp), DIMENSION(SIZE(dbcs)+1) :: alpha, alpha_rotd, xlb, xub, ylb, yub, &
818 24 : zlb, zub
819 :
820 12 : CALL timeset(routineN, handle)
821 :
822 12 : logger => cp_get_default_logger()
823 12 : IF (logger%para_env%is_source()) THEN
824 6 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
825 : ELSE
826 : unit_nr = -1
827 : END IF
828 :
829 12 : Lx = cell_xtnts(2, 1) - cell_xtnts(1, 1)
830 12 : Ly = cell_xtnts(2, 2) - cell_xtnts(1, 2)
831 12 : Lz = cell_xtnts(2, 3) - cell_xtnts(1, 3)
832 :
833 20 : SELECT CASE (parallel_axis)
834 : CASE (x_axis)
835 8 : IF (xtnt(1) < cell_xtnts(1, 1) .OR. xtnt(2) > cell_xtnts(2, 1)) THEN
836 : CALL cp_abort(__LOCATION__, &
837 : "The length of the cylindrical Dirichlet region is larger than the "// &
838 0 : "x range of the simulation cell.")
839 : END IF
840 8 : forb_xtnt1 = base_center(1) - base_radius < cell_xtnts(1, 2)
841 8 : forb_xtnt2 = base_center(1) + base_radius > cell_xtnts(2, 2)
842 8 : forb_xtnt3 = base_center(2) - base_radius < cell_xtnts(1, 3)
843 8 : forb_xtnt4 = base_center(2) + base_radius > cell_xtnts(2, 3)
844 8 : IF (forb_xtnt1 .OR. forb_xtnt2 .OR. forb_xtnt3 .OR. forb_xtnt4) THEN
845 0 : CPABORT("The cylinder does not fit entirely inside the simulation cell.")
846 : END IF
847 : CASE (y_axis)
848 2 : IF (xtnt(1) < cell_xtnts(1, 2) .OR. xtnt(2) > cell_xtnts(2, 2)) THEN
849 : CALL cp_abort(__LOCATION__, &
850 : "The length of the cylindrical Dirichlet region is larger than the "// &
851 0 : "y range of the simulation cell.")
852 : END IF
853 2 : forb_xtnt1 = base_center(1) - base_radius < cell_xtnts(1, 1)
854 2 : forb_xtnt2 = base_center(1) + base_radius > cell_xtnts(2, 1)
855 2 : forb_xtnt3 = base_center(2) - base_radius < cell_xtnts(1, 3)
856 2 : forb_xtnt4 = base_center(2) + base_radius > cell_xtnts(2, 3)
857 2 : IF (forb_xtnt1 .OR. forb_xtnt2 .OR. forb_xtnt3 .OR. forb_xtnt4) THEN
858 0 : CPABORT("The cylinder does not fit entirely inside the simulation cell.")
859 : END IF
860 : CASE (z_axis)
861 2 : IF (xtnt(1) < cell_xtnts(1, 3) .OR. xtnt(2) > cell_xtnts(2, 3)) THEN
862 : CALL cp_abort(__LOCATION__, &
863 : "The length of the cylindrical Dirichlet region is larger than the "// &
864 0 : "z range of the simulation cell.")
865 : END IF
866 2 : forb_xtnt1 = base_center(1) - base_radius < cell_xtnts(1, 1)
867 2 : forb_xtnt2 = base_center(1) + base_radius > cell_xtnts(2, 1)
868 2 : forb_xtnt3 = base_center(2) - base_radius < cell_xtnts(1, 2)
869 2 : forb_xtnt4 = base_center(2) + base_radius > cell_xtnts(2, 2)
870 14 : IF (forb_xtnt1 .OR. forb_xtnt2 .OR. forb_xtnt3 .OR. forb_xtnt4) THEN
871 0 : CPABORT("The cylinder does not fit entirely inside the simulation cell.")
872 : END IF
873 : END SELECT
874 :
875 12 : intern_apx_type = CIRCUMSCRIBED
876 12 : IF (PRESENT(apx_type)) intern_apx_type = apx_type
877 :
878 12 : n_dbcs = SIZE(dbcs)
879 :
880 12 : cntrlaxis_lb = xtnt(1)
881 12 : cntrlaxis_ub = xtnt(2)
882 :
883 12 : theta = twopi/n_dbcs
884 :
885 12 : IF (intern_apx_type == INSCRIBED) THEN
886 6 : h = base_radius ! inscribed uniform prism
887 6 : ELSE IF (intern_apx_type == CIRCUMSCRIBED) THEN
888 6 : h = base_radius/COS(0.5*theta) ! circumscribed uniform prism
889 : ELSE
890 0 : CPABORT("Unknown approximation type for cylinder.")
891 : END IF
892 8 : SELECT CASE (parallel_axis)
893 : CASE (x_axis)
894 32 : IF (h > MINVAL([Ly, Lz])) THEN
895 0 : CPABORT("Reduce the base radius!")
896 : END IF
897 : CASE (y_axis)
898 8 : IF (h > MINVAL([Lx, Lz])) THEN
899 0 : CPABORT("Reduce the base radius!")
900 : END IF
901 : CASE (z_axis)
902 18 : IF (h > MINVAL([Lx, Ly])) THEN
903 0 : CPABORT("Reduce the base radius!")
904 : END IF
905 : END SELECT
906 :
907 12 : alpha = linspace(0.0_dp, 2*pi, n_dbcs + 1)
908 116 : alpha_rotd = alpha + delta_alpha;
909 8 : SELECT CASE (parallel_axis)
910 : CASE (x_axis)
911 72 : DO j = 1, n_dbcs
912 64 : ylb(j) = base_center(1) + h*SIN(alpha(j))
913 64 : zlb(j) = base_center(2) + h*COS(alpha(j))
914 64 : yub(j) = base_center(1) + h*SIN(alpha_rotd(j))
915 72 : zub(j) = base_center(2) + h*COS(alpha_rotd(j))
916 : END DO
917 8 : ylb(n_dbcs + 1) = ylb(1)
918 8 : yub(n_dbcs + 1) = yub(1)
919 8 : zlb(n_dbcs + 1) = zlb(1)
920 8 : zub(n_dbcs + 1) = zub(1)
921 72 : DO j = 1, n_dbcs
922 2112 : ALLOCATE (dbcs(j)%dirichlet_bc)
923 64 : dbcs(j)%dirichlet_bc%dbc_geom = CYLINDRICAL
924 64 : dbcs(j)%dirichlet_bc%v_D = v_D
925 64 : dbcs(j)%dirichlet_bc%osc_frac = osc_frac
926 64 : dbcs(j)%dirichlet_bc%frequency = frequency
927 64 : dbcs(j)%dirichlet_bc%phase = phase
928 64 : dbcs(j)%dirichlet_bc%dbc_id = CYLINDRICAL + j
929 64 : dbcs(j)%dirichlet_bc%smoothing_width = sigma
930 :
931 256 : A = [cntrlaxis_lb, yub(j), zub(j)]
932 256 : B = [cntrlaxis_lb, ylb(j + 1), zlb(j + 1)]
933 256 : C = [cntrlaxis_ub, ylb(j + 1), zlb(j + 1)]
934 256 : D = [cntrlaxis_ub, yub(j), zub(j)]
935 448 : normal_vec = vector_product((A - C), (D - B))
936 448 : unit_normal = normal_vec/SQRT(SUM(normal_vec**2))
937 :
938 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 1) = A
939 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 2) = B
940 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 3) = C
941 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 4) = D
942 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 5) = D + thickness*unit_normal
943 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 6) = A + thickness*unit_normal
944 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 7) = B + thickness*unit_normal
945 256 : dbcs(j)%dirichlet_bc%vertices(1:3, 8) = C + thickness*unit_normal
946 :
947 64 : dbcs(j)%dirichlet_bc%n_tiles = 1
948 :
949 64 : IF ((unit_nr > 0) .AND. verbose) THEN
950 0 : WRITE (unit_nr, '(T3,A,A)') "======== verbose ", REPEAT('=', 61)
951 0 : WRITE (unit_nr, '(T7,A,I5,T20,A,I5,A)') "edge", j, "of the", n_dbcs, &
952 0 : "-gonal prism approximating the cylinder"
953 0 : DO i = 1, 8
954 0 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, &
955 0 : angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
956 : END DO
957 : END IF
958 :
959 : CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
960 72 : x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
961 : END DO
962 : CASE (y_axis)
963 16 : DO j = 1, n_dbcs
964 14 : xlb(j) = base_center(1) + h*SIN(alpha(j))
965 14 : zlb(j) = base_center(2) + h*COS(alpha(j))
966 14 : xub(j) = base_center(1) + h*SIN(alpha_rotd(j))
967 16 : zub(j) = base_center(2) + h*COS(alpha_rotd(j))
968 : END DO
969 2 : xlb(n_dbcs + 1) = xlb(1)
970 2 : xub(n_dbcs + 1) = xub(1)
971 2 : zlb(n_dbcs + 1) = zlb(1)
972 2 : zub(n_dbcs + 1) = zub(1)
973 16 : DO j = 1, n_dbcs
974 462 : ALLOCATE (dbcs(j)%dirichlet_bc)
975 14 : dbcs(j)%dirichlet_bc%dbc_geom = CYLINDRICAL
976 14 : dbcs(j)%dirichlet_bc%v_D = v_D
977 14 : dbcs(j)%dirichlet_bc%osc_frac = osc_frac
978 14 : dbcs(j)%dirichlet_bc%frequency = frequency
979 14 : dbcs(j)%dirichlet_bc%phase = phase
980 14 : dbcs(j)%dirichlet_bc%dbc_id = CYLINDRICAL + j
981 14 : dbcs(j)%dirichlet_bc%smoothing_width = sigma
982 :
983 56 : A = [xub(j), cntrlaxis_lb, zub(j)]
984 56 : B = [xlb(j + 1), cntrlaxis_lb, zlb(j + 1)]
985 56 : C = [xlb(j + 1), cntrlaxis_ub, zlb(j + 1)]
986 56 : D = [xub(j), cntrlaxis_ub, zub(j)]
987 98 : normal_vec = vector_product((A - C), (D - B))
988 98 : unit_normal = normal_vec/SQRT(SUM(normal_vec**2))
989 :
990 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 1) = A
991 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 2) = B
992 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 3) = C
993 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 4) = D
994 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 5) = D + thickness*unit_normal
995 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 6) = A + thickness*unit_normal
996 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 7) = B + thickness*unit_normal
997 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 8) = C + thickness*unit_normal
998 :
999 14 : dbcs(j)%dirichlet_bc%n_tiles = 1
1000 :
1001 14 : IF ((unit_nr > 0) .AND. verbose) THEN
1002 0 : WRITE (unit_nr, '(T3,A,A)') "======== verbose ", REPEAT('=', 61)
1003 0 : WRITE (unit_nr, '(T7,A,I5,T20,A,I5,A)') "edge", j, "of the", n_dbcs, &
1004 0 : "-gonal prism approximating the cylinder"
1005 0 : DO i = 1, 8
1006 0 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, &
1007 0 : angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
1008 : END DO
1009 : END IF
1010 :
1011 : CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
1012 16 : x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
1013 : END DO
1014 : CASE (z_axis)
1015 16 : DO j = 1, n_dbcs
1016 14 : xlb(j) = base_center(1) + h*SIN(alpha(j))
1017 14 : ylb(j) = base_center(2) + h*COS(alpha(j))
1018 14 : xub(j) = base_center(1) + h*SIN(alpha_rotd(j))
1019 16 : yub(j) = base_center(2) + h*COS(alpha_rotd(j))
1020 : END DO
1021 2 : xlb(n_dbcs + 1) = xlb(1)
1022 2 : xub(n_dbcs + 1) = xub(1)
1023 2 : ylb(n_dbcs + 1) = ylb(1)
1024 2 : yub(n_dbcs + 1) = yub(1)
1025 28 : DO j = 1, n_dbcs
1026 462 : ALLOCATE (dbcs(j)%dirichlet_bc)
1027 14 : dbcs(j)%dirichlet_bc%dbc_geom = CYLINDRICAL
1028 14 : dbcs(j)%dirichlet_bc%v_D = v_D
1029 14 : dbcs(j)%dirichlet_bc%osc_frac = osc_frac
1030 14 : dbcs(j)%dirichlet_bc%frequency = frequency
1031 14 : dbcs(j)%dirichlet_bc%phase = phase
1032 14 : dbcs(j)%dirichlet_bc%dbc_id = CYLINDRICAL + j
1033 14 : dbcs(j)%dirichlet_bc%smoothing_width = sigma
1034 :
1035 56 : A = [xub(j), yub(j), cntrlaxis_lb]
1036 56 : B = [xlb(j + 1), ylb(j + 1), cntrlaxis_lb]
1037 56 : C = [xlb(j + 1), ylb(j + 1), cntrlaxis_ub]
1038 56 : D = [xub(j), yub(j), cntrlaxis_ub]
1039 98 : normal_vec = vector_product((A - C), (D - B))
1040 98 : unit_normal = normal_vec/SQRT(SUM(normal_vec**2))
1041 :
1042 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 1) = A
1043 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 2) = B
1044 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 3) = C
1045 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 4) = D
1046 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 5) = D + thickness*unit_normal
1047 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 6) = A + thickness*unit_normal
1048 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 7) = B + thickness*unit_normal
1049 56 : dbcs(j)%dirichlet_bc%vertices(1:3, 8) = C + thickness*unit_normal
1050 :
1051 14 : dbcs(j)%dirichlet_bc%n_tiles = 1
1052 :
1053 14 : IF ((unit_nr > 0) .AND. verbose) THEN
1054 0 : WRITE (unit_nr, '(T3,A,A)') "======== verbose ", REPEAT('=', 61)
1055 0 : WRITE (unit_nr, '(T7,A,I5,T20,A,I5,A)') "edge", j, "of the", n_dbcs, &
1056 0 : "-gonal prism approximating the cylinder"
1057 0 : DO i = 1, 8
1058 0 : WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, &
1059 0 : angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
1060 : END DO
1061 : END IF
1062 :
1063 : CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
1064 16 : x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
1065 : END DO
1066 : END SELECT
1067 :
1068 12 : CALL timestop(handle)
1069 :
1070 12 : END SUBROUTINE aa_cylindrical_dbc_setup
1071 :
1072 : ! **************************************************************************************************
1073 : !> \brief Creteas an evenly-spaced 1D array between two values
1074 : !> \param a starting value
1075 : !> \param b end value
1076 : !> \param N number of evenly-spaced points between a and b
1077 : !> \return array containg N evenly-spaced points between a and b
1078 : !> \par History
1079 : !> 07.2014 created [Hossein Bani-Hashemian]
1080 : !> \author Mohammad Hossein Bani-Hashemian
1081 : ! **************************************************************************************************
1082 12 : FUNCTION linspace(a, b, N) RESULT(arr)
1083 :
1084 : REAL(dp), INTENT(IN) :: a, b
1085 : INTEGER, INTENT(IN) :: N
1086 : REAL(dp), DIMENSION(N) :: arr
1087 :
1088 : INTEGER :: i
1089 : REAL(dp) :: dx
1090 :
1091 12 : dx = (b - a)/(N - 1)
1092 232 : arr = [(a + (i - 1)*dx, i=1, N)]
1093 :
1094 12 : END FUNCTION linspace
1095 :
1096 : ! **************************************************************************************************
1097 : !> \brief rotates and translates a cuboid
1098 : !> \param cuboid_vtx vertices of the cuboid
1099 : !> \param Rmat rotation matrix
1100 : !> \param Tpnt translation vector (translates the origin to the center of the cuboid)
1101 : !> \param cuboid_transfd_vtx vertices of the transformed cuboid
1102 : !> \par History
1103 : !> 11.2015 created [Hossein Bani-Hashemian]
1104 : !> \author Mohammad Hossein Bani-Hashemian
1105 : ! **************************************************************************************************
1106 160 : SUBROUTINE rotate_translate_cuboid(cuboid_vtx, Rmat, Tpnt, cuboid_transfd_vtx)
1107 :
1108 : REAL(dp), DIMENSION(3, 8), INTENT(IN) :: cuboid_vtx
1109 : REAL(dp), DIMENSION(3, 3), INTENT(OUT) :: Rmat
1110 : REAL(dp), DIMENSION(3), INTENT(OUT) :: Tpnt
1111 : REAL(dp), DIMENSION(3, 8), INTENT(OUT) :: cuboid_transfd_vtx
1112 :
1113 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rotate_translate_cuboid'
1114 :
1115 : INTEGER :: handle
1116 : REAL(dp), DIMENSION(3) :: A, A2, AA2, AB, AD, B, B2, C, C2, D, D2, &
1117 : e1_locl, e2_locl, e3_locl
1118 :
1119 160 : CALL timeset(routineN, handle)
1120 :
1121 160 : Rmat = 0.0_dp
1122 :
1123 1120 : A = cuboid_vtx(1:3, 1); A2 = cuboid_vtx(1:3, 6)
1124 1120 : B = cuboid_vtx(1:3, 2); B2 = cuboid_vtx(1:3, 7)
1125 1120 : C = cuboid_vtx(1:3, 3); C2 = cuboid_vtx(1:3, 8)
1126 1120 : D = cuboid_vtx(1:3, 4); D2 = cuboid_vtx(1:3, 5)
1127 :
1128 640 : Tpnt = (A + C2)/2.0_dp
1129 :
1130 640 : AB = B - A
1131 640 : AD = D - A
1132 640 : AA2 = A2 - A
1133 :
1134 : ! unit vectors generating the local coordinate system
1135 1120 : e1_locl = AB/SQRT(SUM(AB**2))
1136 1120 : e2_locl = AD/SQRT(SUM(AD**2))
1137 1120 : e3_locl = AA2/SQRT(SUM(AA2**2))
1138 : ! rotation matrix
1139 640 : Rmat(1:3, 1) = e1_locl
1140 640 : Rmat(1:3, 2) = e2_locl
1141 640 : Rmat(1:3, 3) = e3_locl
1142 :
1143 160 : cuboid_transfd_vtx(1:3, 1) = rotate_translate_vector(Rmat, Tpnt, BWROT, A)
1144 160 : cuboid_transfd_vtx(1:3, 2) = rotate_translate_vector(Rmat, Tpnt, BWROT, B)
1145 160 : cuboid_transfd_vtx(1:3, 3) = rotate_translate_vector(Rmat, Tpnt, BWROT, C)
1146 160 : cuboid_transfd_vtx(1:3, 4) = rotate_translate_vector(Rmat, Tpnt, BWROT, D)
1147 160 : cuboid_transfd_vtx(1:3, 5) = rotate_translate_vector(Rmat, Tpnt, BWROT, D2)
1148 160 : cuboid_transfd_vtx(1:3, 6) = rotate_translate_vector(Rmat, Tpnt, BWROT, A2)
1149 160 : cuboid_transfd_vtx(1:3, 7) = rotate_translate_vector(Rmat, Tpnt, BWROT, B2)
1150 160 : cuboid_transfd_vtx(1:3, 8) = rotate_translate_vector(Rmat, Tpnt, BWROT, C2)
1151 :
1152 160 : CALL timestop(handle)
1153 :
1154 160 : END SUBROUTINE rotate_translate_cuboid
1155 :
1156 : ! **************************************************************************************************
1157 : !> \brief transforms a position vector according to a given rotation matrix and a
1158 : !> translation vector
1159 : !> \param Rmat rotation matrix
1160 : !> \param Tp image of the origin under the translation
1161 : !> \param direction forward or backward transformation
1162 : !> \param vec input vector
1163 : !> \return transformed vector
1164 : !> \par History
1165 : !> 11.2015 created [Hossein Bani-Hashemian]
1166 : !> \author Mohammad Hossein Bani-Hashemian
1167 : ! **************************************************************************************************
1168 15291164 : FUNCTION rotate_translate_vector(Rmat, Tp, direction, vec) RESULT(vec_transfd)
1169 : REAL(dp), DIMENSION(3, 3), INTENT(IN) :: Rmat
1170 : REAL(dp), DIMENSION(3), INTENT(IN) :: Tp
1171 : INTEGER, INTENT(IN) :: direction
1172 : REAL(dp), DIMENSION(3), INTENT(IN) :: vec
1173 : REAL(dp), DIMENSION(3) :: vec_transfd
1174 :
1175 : REAL(dp), DIMENSION(3) :: Tpoint, vec_tmp
1176 : REAL(dp), DIMENSION(3, 3) :: Rmat_inv
1177 :
1178 15291164 : Tpoint = 0.0_dp
1179 :
1180 15291164 : IF (direction == FWROT) THEN
1181 0 : vec_tmp = MATMUL(Rmat, vec)
1182 0 : vec_transfd = vec_tmp + Tp
1183 15291164 : ELSEIF (direction == BWROT) THEN
1184 15291164 : Rmat_inv = inv_3x3(Rmat)
1185 15291164 : Tpoint(1) = Rmat_inv(1, 1)*Tp(1) + Rmat_inv(1, 2)*Tp(2) + Rmat_inv(1, 3)*Tp(3)
1186 15291164 : Tpoint(2) = Rmat_inv(2, 1)*Tp(1) + Rmat_inv(2, 2)*Tp(2) + Rmat_inv(2, 3)*Tp(3)
1187 15291164 : Tpoint(3) = Rmat_inv(3, 1)*Tp(1) + Rmat_inv(3, 2)*Tp(2) + Rmat_inv(3, 3)*Tp(3)
1188 198785132 : vec_tmp = MATMUL(Rmat_inv, vec)
1189 61164656 : vec_transfd = vec_tmp - Tpoint
1190 : END IF
1191 :
1192 15291164 : END FUNCTION rotate_translate_vector
1193 :
1194 : ! **************************************************************************************************
1195 : !> \brief Partitions an axis-aligned cuboid into tiles.
1196 : !> \param dbc_vertices vertices of the axis-aligned Dirichlet region to be partitioned
1197 : !> \param n_prtn vector of size 3 specifying the number of times that the x, y and
1198 : !> z interval (defining the region) should be partitioned into
1199 : !> \param tiles the obtained tiles
1200 : !> \par History
1201 : !> 10.2015 created [Hossein Bani-Hashemian]
1202 : !> \author Mohammad Hossein Bani-Hashemian
1203 : ! **************************************************************************************************
1204 30 : SUBROUTINE aa_dbc_partition(dbc_vertices, n_prtn, tiles)
1205 :
1206 : REAL(dp), DIMENSION(3, 8), INTENT(IN) :: dbc_vertices
1207 : INTEGER, DIMENSION(3), INTENT(IN) :: n_prtn
1208 : TYPE(tile_p_type), DIMENSION(:), INTENT(OUT), &
1209 : POINTER :: tiles
1210 :
1211 : CHARACTER(LEN=*), PARAMETER :: routineN = 'aa_dbc_partition'
1212 :
1213 : INTEGER :: handle, i, ii, jj, k, kk
1214 : REAL(dp) :: step
1215 30 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xprtn_lb, xprtn_ub, yprtn_lb, yprtn_ub, &
1216 30 : zprtn_lb, zprtn_ub
1217 : REAL(dp), DIMENSION(2) :: x_xtnt, y_xtnt, z_xtnt
1218 :
1219 30 : CALL timeset(routineN, handle)
1220 :
1221 : ! vertices are labeled as:
1222 : ! 6-------5 A2------D2
1223 : ! /| /| /| /|
1224 : ! 7-|-----8 | B2-|----C2 |
1225 : ! | 1-----|-4 or | A-----|-D
1226 : ! |/ |/ |/ |/
1227 : ! 2-------3 B-------C
1228 :
1229 120 : ALLOCATE (xprtn_lb(n_prtn(1)), xprtn_ub(n_prtn(1)))
1230 120 : ALLOCATE (yprtn_lb(n_prtn(2)), yprtn_ub(n_prtn(2)))
1231 120 : ALLOCATE (zprtn_lb(n_prtn(3)), zprtn_ub(n_prtn(3)))
1232 :
1233 : ! find the extents of the cuboid in all three directions
1234 600 : x_xtnt(1) = MINVAL(dbc_vertices(1, :)); x_xtnt(2) = MAXVAL(dbc_vertices(1, :))
1235 600 : y_xtnt(1) = MINVAL(dbc_vertices(2, :)); y_xtnt(2) = MAXVAL(dbc_vertices(2, :))
1236 600 : z_xtnt(1) = MINVAL(dbc_vertices(3, :)); z_xtnt(2) = MAXVAL(dbc_vertices(3, :))
1237 :
1238 : ! divide the x, y and z extents into n_prtn partitions
1239 30 : step = (x_xtnt(2) - x_xtnt(1))/REAL(n_prtn(1), kind=dp)
1240 140 : xprtn_lb(:) = x_xtnt(1) + [(i, i=0, n_prtn(1) - 1)]*step ! lower bounds
1241 140 : xprtn_ub(:) = x_xtnt(1) + [(i, i=1, n_prtn(1))]*step ! upper bounds
1242 :
1243 30 : step = (y_xtnt(2) - y_xtnt(1))/REAL(n_prtn(2), kind=dp)
1244 136 : yprtn_lb(:) = y_xtnt(1) + [(i, i=0, n_prtn(2) - 1)]*step
1245 136 : yprtn_ub(:) = y_xtnt(1) + [(i, i=1, n_prtn(2))]*step
1246 :
1247 30 : step = (z_xtnt(2) - z_xtnt(1))/REAL(n_prtn(3), kind=dp)
1248 140 : zprtn_lb(:) = z_xtnt(1) + [(i, i=0, n_prtn(3) - 1)]*step
1249 140 : zprtn_ub(:) = z_xtnt(1) + [(i, i=1, n_prtn(3))]*step
1250 :
1251 154 : ALLOCATE (tiles(n_prtn(1)*n_prtn(2)*n_prtn(3)))
1252 : k = 1
1253 70 : DO kk = 1, n_prtn(3)
1254 118 : DO jj = 1, n_prtn(2)
1255 152 : DO ii = 1, n_prtn(1)
1256 2176 : ALLOCATE (tiles(k)%tile)
1257 64 : tiles(k)%tile%tile_id = k
1258 :
1259 256 : tiles(k)%tile%vertices(1:3, 1) = [xprtn_lb(ii), yprtn_lb(jj), zprtn_lb(kk)]
1260 256 : tiles(k)%tile%vertices(1:3, 2) = [xprtn_ub(ii), yprtn_lb(jj), zprtn_lb(kk)]
1261 256 : tiles(k)%tile%vertices(1:3, 3) = [xprtn_ub(ii), yprtn_ub(jj), zprtn_lb(kk)]
1262 256 : tiles(k)%tile%vertices(1:3, 4) = [xprtn_lb(ii), yprtn_ub(jj), zprtn_lb(kk)]
1263 256 : tiles(k)%tile%vertices(1:3, 5) = [xprtn_lb(ii), yprtn_ub(jj), zprtn_ub(kk)]
1264 256 : tiles(k)%tile%vertices(1:3, 6) = [xprtn_lb(ii), yprtn_lb(jj), zprtn_ub(kk)]
1265 256 : tiles(k)%tile%vertices(1:3, 7) = [xprtn_ub(ii), yprtn_lb(jj), zprtn_ub(kk)]
1266 256 : tiles(k)%tile%vertices(1:3, 8) = [xprtn_ub(ii), yprtn_ub(jj), zprtn_ub(kk)]
1267 :
1268 64 : tiles(k)%tile%volume = 0.0_dp
1269 112 : k = k + 1
1270 : END DO
1271 : END DO
1272 : END DO
1273 :
1274 30 : CALL timestop(handle)
1275 :
1276 60 : END SUBROUTINE aa_dbc_partition
1277 :
1278 : ! **************************************************************************************************
1279 : !> \brief Partitions an arbitrary cuboid into tiles.
1280 : !> \param dbc_vertices vertices of the axis-aligned Dirichlet region to be partitioned
1281 : !> \param n_prtn number of partitions along the edges
1282 : !> \param tiles the obtained tiles
1283 : !> \par History
1284 : !> 10.2015 created [Hossein Bani-Hashemian]
1285 : !> \author Mohammad Hossein Bani-Hashemian
1286 : ! **************************************************************************************************
1287 98 : SUBROUTINE arbitrary_dbc_partition(dbc_vertices, n_prtn, tiles)
1288 :
1289 : REAL(dp), DIMENSION(3, 8), INTENT(IN) :: dbc_vertices
1290 : INTEGER, DIMENSION(2), INTENT(IN) :: n_prtn
1291 : TYPE(tile_p_type), DIMENSION(:), INTENT(OUT), &
1292 : POINTER :: tiles
1293 :
1294 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arbitrary_dbc_partition'
1295 :
1296 : INTEGER :: handle, i, k, l, np1, np2
1297 : REAL(dp) :: ABlength, ADlength, step, thickness
1298 98 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: end_points, X
1299 : REAL(dp), DIMENSION(3) :: A, AB, AC, AD, B, C, D, D2, &
1300 : normal_vector, point1, point2, &
1301 : unit_normal
1302 :
1303 98 : CALL timeset(routineN, handle)
1304 :
1305 392 : A = dbc_vertices(:, 1)
1306 392 : B = dbc_vertices(:, 2)
1307 392 : C = dbc_vertices(:, 3)
1308 392 : D = dbc_vertices(:, 4)
1309 392 : D2 = dbc_vertices(:, 5)
1310 :
1311 392 : AB = B - A
1312 392 : AC = C - A
1313 392 : AD = D - A
1314 98 : normal_vector = vector_product(AB, AC)
1315 686 : unit_normal = normal_vector/SQRT(SUM(normal_vector**2))
1316 392 : thickness = SQRT(SUM((D - D2)**2))
1317 :
1318 : ! the larger n_prtn is assigned to the longer edge
1319 392 : ABlength = SQRT(SUM(AB**2))
1320 392 : ADlength = SQRT(SUM(AD**2))
1321 98 : IF (ADlength >= ABlength) THEN
1322 92 : np1 = MAX(n_prtn(1) + 1, n_prtn(2) + 1)
1323 92 : np2 = MIN(n_prtn(1) + 1, n_prtn(2) + 1)
1324 : ELSE
1325 6 : np1 = MIN(n_prtn(1) + 1, n_prtn(2) + 1)
1326 6 : np2 = MAX(n_prtn(1) + 1, n_prtn(2) + 1)
1327 : END IF
1328 :
1329 490 : ALLOCATE (X(np1, np2, 3))
1330 392 : ALLOCATE (end_points(2, np1, 3))
1331 :
1332 : ! partition AD and BC
1333 346 : DO l = 1, np1
1334 248 : step = REAL((l - 1), kind=dp)/REAL((np1 - 1), kind=dp)
1335 992 : end_points(1, l, :) = A*(1.0_dp - step) + D*step
1336 1090 : end_points(2, l, :) = B*(1.0_dp - step) + C*step
1337 : END DO
1338 :
1339 : ! partition in the second direction along the line segments with endpoints from
1340 : ! the previous partitioning step
1341 346 : DO l = 1, np1
1342 992 : point1(:) = end_points(1, l, :)
1343 992 : point2(:) = end_points(2, l, :)
1344 858 : DO i = 1, np2
1345 512 : step = REAL((i - 1), kind=dp)/REAL((np2 - 1), kind=dp)
1346 2296 : X(l, i, :) = point1*(1.0_dp - step) + point2*step
1347 : END DO
1348 : END DO
1349 :
1350 454 : ALLOCATE (tiles((np1 - 1)*(np2 - 1)))
1351 : k = 1
1352 248 : DO l = 1, np1 - 1
1353 408 : DO i = 1, np2 - 1
1354 5280 : ALLOCATE (tiles(k)%tile)
1355 160 : tiles(k)%tile%tile_id = k
1356 :
1357 640 : tiles(k)%tile%vertices(1:3, 1) = X(l, i, :)
1358 640 : tiles(k)%tile%vertices(1:3, 2) = X(l, i + 1, :)
1359 640 : tiles(k)%tile%vertices(1:3, 3) = X(l + 1, i + 1, :)
1360 640 : tiles(k)%tile%vertices(1:3, 4) = X(l + 1, i, :)
1361 640 : tiles(k)%tile%vertices(1:3, 5) = X(l + 1, i, :) + thickness*unit_normal
1362 640 : tiles(k)%tile%vertices(1:3, 6) = X(l, i, :) + thickness*unit_normal
1363 640 : tiles(k)%tile%vertices(1:3, 7) = X(l, i + 1, :) + thickness*unit_normal
1364 640 : tiles(k)%tile%vertices(1:3, 8) = X(l + 1, i + 1, :) + thickness*unit_normal
1365 :
1366 160 : tiles(k)%tile%volume = 0.0_dp
1367 310 : k = k + 1
1368 : END DO
1369 : END DO
1370 :
1371 98 : CALL timestop(handle)
1372 :
1373 196 : END SUBROUTINE arbitrary_dbc_partition
1374 :
1375 : ! **************************************************************************************************
1376 : !> \brief computes the pw data corresponding to an axis-aligned tile
1377 : !> \param cell_xtnts extents of the simulation cell
1378 : !> \param x_locl x grid vector of the simulation box local to this process
1379 : !> \param y_locl y grid vector of the simulation box local to this process
1380 : !> \param z_locl z grid vector of the simulation box local to this process
1381 : !> \param tile_vertices coordinates of the vertices of the input tile
1382 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
1383 : !> \param is_periodic whether or not the tile is periodic with a periodicity of one unit cell
1384 : !> \param tile_pw the output pw data
1385 : !> \par History
1386 : !> 10.2015 created [Hossein Bani-Hashemian]
1387 : !> \author Mohammad Hossein Bani-Hashemian
1388 : ! **************************************************************************************************
1389 64 : SUBROUTINE aa_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, tile_vertices, &
1390 : sigma, is_periodic, tile_pw)
1391 :
1392 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
1393 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
1394 : REAL(dp), DIMENSION(3, 8), INTENT(IN) :: tile_vertices
1395 : REAL(dp), INTENT(IN) :: sigma
1396 : LOGICAL, INTENT(IN) :: is_periodic
1397 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: tile_pw
1398 :
1399 : CHARACTER(LEN=*), PARAMETER :: routineN = 'aa_tile_pw_compute'
1400 :
1401 : INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
1402 : ub2, ub3
1403 : INTEGER, DIMENSION(2, 3) :: bounds_local
1404 : REAL(dp) :: fx, fy, fz, gx, gy, gz, hx, hy, hz, xi0, &
1405 : xil, xir, yj0, yjl, yjr, zk0, zkl, zkr
1406 64 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xlocll, xloclr, ylocll, yloclr, zlocll, &
1407 64 : zloclr
1408 : REAL(dp), DIMENSION(2) :: x_xtnt, y_xtnt, z_xtnt
1409 : REAL(dp), DIMENSION(3) :: per
1410 : REAL(dp), DIMENSION(3, 3) :: prefactor
1411 :
1412 64 : CALL timeset(routineN, handle)
1413 :
1414 640 : bounds_local = tile_pw%pw_grid%bounds_local
1415 64 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1416 64 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1417 64 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1418 :
1419 448 : ALLOCATE (xlocll(lb1:ub1), ylocll(lb2:ub2), zlocll(lb3:ub3))
1420 256 : ALLOCATE (xloclr(lb1:ub1), yloclr(lb2:ub2), zloclr(lb3:ub3))
1421 :
1422 : ! periodicities
1423 256 : per = cell_xtnts(2, :) - cell_xtnts(1, :)
1424 :
1425 1216 : x_xtnt(1) = MINVAL(tile_vertices(1, :)); x_xtnt(2) = MAXVAL(tile_vertices(1, :))
1426 1216 : y_xtnt(1) = MINVAL(tile_vertices(2, :)); y_xtnt(2) = MAXVAL(tile_vertices(2, :))
1427 1216 : z_xtnt(1) = MINVAL(tile_vertices(3, :)); z_xtnt(2) = MAXVAL(tile_vertices(3, :))
1428 :
1429 : ! prefactor: columns: left, original, right - rows: x , y, z directions
1430 832 : prefactor = 0.5_dp
1431 :
1432 64 : IF (is_periodic) THEN
1433 0 : DO k = lb3, ub3
1434 0 : DO j = lb2, ub2
1435 0 : DO i = lb1, ub1
1436 0 : xi0 = x_locl(i); yj0 = y_locl(j); zk0 = z_locl(k)
1437 0 : xil = x_locl(i) - per(1); yjl = y_locl(j) - per(2); zkl = z_locl(k) - per(3)
1438 0 : xir = x_locl(i) + per(1); yjr = y_locl(j) + per(2); zkr = z_locl(k) + per(3)
1439 :
1440 : ! points from the original cell local to the current processor
1441 0 : fx = prefactor(1, 2)*(erf((xi0 - x_xtnt(1))/sigma) - erf((xi0 - x_xtnt(2))/sigma))
1442 0 : fy = prefactor(2, 2)*(erf((yj0 - y_xtnt(1))/sigma) - erf((yj0 - y_xtnt(2))/sigma))
1443 0 : fz = prefactor(3, 2)*(erf((zk0 - z_xtnt(1))/sigma) - erf((zk0 - z_xtnt(2))/sigma))
1444 :
1445 : ! periodically replicated cell on the left, points local to the current processor
1446 0 : gx = prefactor(1, 1)*(erf((xil - x_xtnt(1))/sigma) - erf((xil - x_xtnt(2))/sigma))
1447 0 : gy = prefactor(2, 1)*(erf((yjl - y_xtnt(1))/sigma) - erf((yjl - y_xtnt(2))/sigma))
1448 0 : gz = prefactor(3, 1)*(erf((zkl - z_xtnt(1))/sigma) - erf((zkl - z_xtnt(2))/sigma))
1449 :
1450 : ! periodically replicated cell on the right, points local to the current processor
1451 0 : hx = prefactor(1, 3)*(erf((xir - x_xtnt(1))/sigma) - erf((xir - x_xtnt(2))/sigma))
1452 0 : hy = prefactor(2, 3)*(erf((yjr - y_xtnt(1))/sigma) - erf((yjr - y_xtnt(2))/sigma))
1453 0 : hz = prefactor(3, 3)*(erf((zkr - z_xtnt(1))/sigma) - erf((zkr - z_xtnt(2))/sigma))
1454 :
1455 0 : tile_pw%array(i, j, k) = (fx + gx + hx)*(fy + gy + hy)*(fz + gz + hz)
1456 : END DO
1457 : END DO
1458 : END DO
1459 : ELSE
1460 3926 : DO k = lb3, ub3
1461 249952 : DO j = lb2, ub2
1462 8269833 : DO i = lb1, ub1
1463 8019945 : xi0 = x_locl(i); yj0 = y_locl(j); zk0 = z_locl(k)
1464 :
1465 8019945 : fx = 0.5_dp*(erf((xi0 - x_xtnt(1))/sigma) - erf((xi0 - x_xtnt(2))/sigma))
1466 8019945 : fy = 0.5_dp*(erf((yj0 - y_xtnt(1))/sigma) - erf((yj0 - y_xtnt(2))/sigma))
1467 8019945 : fz = 0.5_dp*(erf((zk0 - z_xtnt(1))/sigma) - erf((zk0 - z_xtnt(2))/sigma))
1468 :
1469 8265971 : tile_pw%array(i, j, k) = fx*fy*fz
1470 : END DO
1471 : END DO
1472 : END DO
1473 : END IF
1474 :
1475 64 : CALL timestop(handle)
1476 :
1477 128 : END SUBROUTINE aa_tile_pw_compute
1478 :
1479 : ! **************************************************************************************************
1480 : !> \brief computes the pw data corresponding to an arbitrary (rectangular-shaped) tile
1481 : !> \param cell_xtnts extents of the simulation cell
1482 : !> \param x_locl x grid vector of the simulation box local to this process
1483 : !> \param y_locl y grid vector of the simulation box local to this process
1484 : !> \param z_locl z grid vector of the simulation box local to this process
1485 : !> \param tile_vertices coordinates of the vertices of the input tile
1486 : !> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
1487 : !> \param tile_pw the output pw data
1488 : !> \par History
1489 : !> 11.2015 created [Hossein Bani-Hashemian]
1490 : !> \author Mohammad Hossein Bani-Hashemian
1491 : ! **************************************************************************************************
1492 160 : SUBROUTINE arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, tile_vertices, &
1493 : sigma, tile_pw)
1494 :
1495 : REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
1496 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
1497 : REAL(dp), DIMENSION(3, 8), INTENT(IN) :: tile_vertices
1498 : REAL(dp), INTENT(IN) :: sigma
1499 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: tile_pw
1500 :
1501 : CHARACTER(LEN=*), PARAMETER :: routineN = 'arbitrary_tile_pw_compute'
1502 :
1503 : INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
1504 : ub2, ub3
1505 : INTEGER, DIMENSION(2, 3) :: bounds_local
1506 : REAL(dp) :: cm_x, cm_y, cm_z, fx, fy, fz, xi0, yj0, &
1507 : zk0
1508 160 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xlocll, xloclr, ylocll, yloclr, zlocll, &
1509 160 : zloclr
1510 : REAL(dp), DIMENSION(2) :: transfdx_xtnt, transfdy_xtnt, &
1511 : transfdz_xtnt, x_xtnt, y_xtnt, z_xtnt
1512 : REAL(dp), DIMENSION(3) :: per, pnt0, Tpnt
1513 : REAL(dp), DIMENSION(3, 3) :: prefactor, Rmat
1514 : REAL(dp), DIMENSION(3, 8) :: tile_transfd_vertices
1515 :
1516 160 : CALL timeset(routineN, handle)
1517 :
1518 1600 : bounds_local = tile_pw%pw_grid%bounds_local
1519 160 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1520 160 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1521 160 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1522 :
1523 1120 : ALLOCATE (xlocll(lb1:ub1), ylocll(lb2:ub2), zlocll(lb3:ub3))
1524 640 : ALLOCATE (xloclr(lb1:ub1), yloclr(lb2:ub2), zloclr(lb3:ub3))
1525 :
1526 640 : per = cell_xtnts(2, :) - cell_xtnts(1, :)
1527 :
1528 160 : CALL rotate_translate_cuboid(tile_vertices, Rmat, Tpnt, tile_transfd_vertices)
1529 :
1530 1600 : transfdx_xtnt(1) = MINVAL(tile_transfd_vertices(1, :))
1531 1600 : transfdx_xtnt(2) = MAXVAL(tile_transfd_vertices(1, :))
1532 1600 : transfdy_xtnt(1) = MINVAL(tile_transfd_vertices(2, :))
1533 1600 : transfdy_xtnt(2) = MAXVAL(tile_transfd_vertices(2, :))
1534 1600 : transfdz_xtnt(1) = MINVAL(tile_transfd_vertices(3, :))
1535 1600 : transfdz_xtnt(2) = MAXVAL(tile_transfd_vertices(3, :))
1536 :
1537 160 : cm_x = 0.5_dp*(transfdx_xtnt(2) + transfdx_xtnt(1))
1538 160 : cm_y = 0.5_dp*(transfdy_xtnt(2) + transfdy_xtnt(1))
1539 160 : cm_z = 0.5_dp*(transfdz_xtnt(2) + transfdz_xtnt(1))
1540 :
1541 480 : x_xtnt = transfdx_xtnt - cm_x
1542 480 : y_xtnt = transfdy_xtnt - cm_y
1543 480 : z_xtnt = transfdz_xtnt - cm_z
1544 :
1545 2080 : prefactor = 0.5_dp
1546 :
1547 9124 : DO k = lb3, ub3
1548 525364 : DO j = lb2, ub2
1549 15815088 : DO i = lb1, ub1
1550 61159536 : pnt0 = rotate_translate_vector(Rmat, Tpnt, BWROT, [x_locl(i), y_locl(j), z_locl(k)])
1551 :
1552 15289884 : xi0 = pnt0(1); yj0 = pnt0(2); zk0 = pnt0(3)
1553 :
1554 15289884 : fx = 0.5_dp*(erf((xi0 - x_xtnt(1))/sigma) - erf((xi0 - x_xtnt(2))/sigma))
1555 15289884 : fy = 0.5_dp*(erf((yj0 - y_xtnt(1))/sigma) - erf((yj0 - y_xtnt(2))/sigma))
1556 15289884 : fz = 0.5_dp*(erf((zk0 - z_xtnt(1))/sigma) - erf((zk0 - z_xtnt(2))/sigma))
1557 :
1558 15806124 : tile_pw%array(i, j, k) = fx*fy*fz
1559 : END DO
1560 : END DO
1561 : END DO
1562 :
1563 160 : CALL timestop(handle)
1564 :
1565 320 : END SUBROUTINE arbitrary_tile_pw_compute
1566 :
1567 : END MODULE dirichlet_bc_methods
|
