Octopus
derivatives.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
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
27
29 use accel_oct_m
31 use batch_oct_m
34 use debug_oct_m
40 use global_oct_m
41 use iso_c_binding
42 use, intrinsic :: iso_fortran_env
45 use loct_oct_m
46 use math_oct_m
47 use mesh_oct_m
53 use parser_oct_m
55 use space_oct_m
63 use types_oct_m
64 use utils_oct_m
66
67! debug purposes
68! use io_binary_oct_m
69! use io_function_oct_m
70! use io_oct_m
71! use unit_oct_m
72! use unit_system_oct_m
73
74 implicit none
75
76 private
77 public :: &
115
116
117 integer, parameter :: &
118 DER_BC_ZERO_F = 0, & !< function is zero at the boundaries
119 der_bc_zero_df = 1, &
120 der_bc_period = 2
122 integer, parameter, public :: &
123 DER_STAR = 1, &
124 der_variational = 2, &
125 der_cube = 3, &
126 der_starplus = 4, &
128
129 integer, parameter :: &
130 BLOCKING = 1, &
131 non_blocking = 2
132
135 type derivatives_t
136 ! Components are public by default
137 type(boundaries_t) :: boundaries
138 type(mesh_t), pointer :: mesh => null()
139 integer :: dim = 0
140 integer, private :: periodic_dim = 0
141 integer :: order = 0
142 integer :: stencil_type = 0
143
144
145 logical, private :: stencil_primitive_coordinates
146
147 real(real64), allocatable :: masses(:)
148
151 real(real64), private :: lapl_cutoff = m_zero
152
153 type(nl_operator_t), allocatable, private :: op(:)
155 type(nl_operator_t), pointer :: lapl => null()
156 type(nl_operator_t), pointer :: grad(:) => null()
157
158 integer, allocatable :: n_ghost(:)
159!#if defined(HAVE_MPI)
160 integer, private :: comm_method = 0
161!#endif
162 type(derivatives_t), pointer :: finer => null()
163 type(derivatives_t), pointer :: coarser => null()
164 type(transfer_table_t), pointer :: to_finer => null()
165 type(transfer_table_t), pointer :: to_coarser => null()
166 end type derivatives_t
167
170 private
171!#ifdef HAVE_MPI
172 type(par_vec_handle_batch_t) :: pv_h
173!#endif
174 type(derivatives_t), pointer :: der
175 type(nl_operator_t), pointer :: op
176 type(batch_t), pointer :: ff
177 type(batch_t), pointer :: opff
178 logical :: ghost_update
179 logical :: factor_present
180 real(real64) :: factor
182
183 type(accel_kernel_t) :: kernel_uvw_xyz, kernel_dcurl, kernel_zcurl
184
185contains
186
187 ! ---------------------------------------------------------
188 subroutine derivatives_init(der, namespace, space, coord_system, order)
189 type(derivatives_t), target, intent(inout) :: der
190 type(namespace_t), intent(in) :: namespace
191 class(space_t), intent(in) :: space
192 class(coordinate_system_t), intent(in) :: coord_system
193 integer, optional, intent(in) :: order
194
195 integer :: idir
196 integer :: default_stencil
197 character(len=40) :: flags
198 logical :: stencil_primitive_coordinates
199
200 push_sub(derivatives_init)
201
202 ! Non-orthogonal curvilinear coordinates are currently not implemented
203 assert(.not. coord_system%local_basis .or. coord_system%orthogonal)
204
205 ! copy this value to my structure
206 der%dim = space%dim
207 der%periodic_dim = space%periodic_dim
208
209 !%Variable DerivativesStencil
210 !%Type integer
211 !%Default stencil_star
212 !%Section Mesh::Derivatives
213 !%Description
214 !% Decides what kind of stencil is used, <i>i.e.</i> which points, around
215 !% each point in the mesh, are the neighboring points used in the
216 !% expression of the differential operator.
217 !%
218 !% If curvilinear coordinates are to be used, then only the <tt>stencil_starplus</tt>
219 !% or the <tt>stencil_cube</tt> may be used. We only recommend the <tt>stencil_starplus</tt>,
220 !% since the cube typically needs far too much memory.
221 !%Option stencil_star 1
222 !% A star around each point (<i>i.e.</i>, only points on the axis).
223 !%Option stencil_variational 2
224 !% Same as the star, but with coefficients built in a different way.
225 !%Option stencil_cube 3
226 !% A cube of points around each point.
227 !%Option stencil_starplus 4
228 !% The star, plus a number of off-axis points.
229 !%Option stencil_stargeneral 5
230 !% The general star. Default for non-orthogonal grids.
231 !%End
232 default_stencil = der_star
233 if (coord_system%local_basis) default_stencil = der_starplus
234 if (.not. coord_system%orthogonal) default_stencil = der_stargeneral
236 call parse_variable(namespace, 'DerivativesStencil', default_stencil, der%stencil_type)
237
238 if (.not. varinfo_valid_option('DerivativesStencil', der%stencil_type)) then
239 call messages_input_error(namespace, 'DerivativesStencil')
240 end if
241 call messages_print_var_option("DerivativesStencil", der%stencil_type, namespace=namespace)
242
243 if (coord_system%local_basis .and. der%stencil_type < der_cube) call messages_input_error(namespace, 'DerivativesStencil')
244 if (der%stencil_type == der_variational) then
245 call parse_variable(namespace, 'DerivativesLaplacianFilter', m_one, der%lapl_cutoff)
246 end if
247
248 !%Variable DerivativesOrder
249 !%Type integer
250 !%Default 4
251 !%Section Mesh::Derivatives
252 !%Description
253 !% This variable gives the discretization order for the approximation of
254 !% the differential operators. This means, basically, that
255 !% <tt>DerivativesOrder</tt> points are used in each positive/negative
256 !% spatial direction, <i>e.g.</i> <tt>DerivativesOrder = 1</tt> would give
257 !% the well-known three-point formula in 1D.
258 !% The number of points actually used for the Laplacian
259 !% depends on the stencil used. Let <math>O</math> = <tt>DerivativesOrder</tt>, and <math>d</math> = <tt>Dimensions</tt>.
260 !% <ul>
261 !% <li> <tt>stencil_star</tt>: <math>2 O d + 1</math>
262 !% <li> <tt>stencil_cube</tt>: <math>(2 O + 1)^d</math>
263 !% <li> <tt>stencil_starplus</tt>: <math>2 O d + 1 + n</math> with <i>n</i> being 8
264 !% in 2D and 24 in 3D.
265 !% </ul>
266 !%End
267 call parse_variable(namespace, 'DerivativesOrder', 4, der%order)
268 ! overwrite order if given as argument
269 if (present(order)) then
270 der%order = order
271 end if
273#ifdef HAVE_MPI
274 !%Variable ParallelizationOfDerivatives
275 !%Type integer
276 !%Default non_blocking
277 !%Section Execution::Parallelization
278 !%Description
279 !% This option selects how the communication of mesh boundaries is performed.
280 !%Option blocking 1
281 !% Blocking communication.
282 !%Option non_blocking 2
283 !% Communication is based on non-blocking point-to-point communication.
284 !%End
285
286 call parse_variable(namespace, 'ParallelizationOfDerivatives', non_blocking, der%comm_method)
287
288 if (.not. varinfo_valid_option('ParallelizationOfDerivatives', der%comm_method)) then
289 call messages_input_error(namespace, 'ParallelizationOfDerivatives')
290 end if
291
292 call messages_obsolete_variable(namespace, 'OverlapDerivatives', 'ParallelizationOfDerivatives')
293#endif
294
295 !%Variable StencilPrimitiveCoordinates
296 !%Type logical
297 !%Section Mesh::Derivatives
298 !%Description
299 !% This variable controls the method for generating the stencil weights for the Laplacian.
300 !% For the gradient, the weights are always computed for primitive coordinates.
301 !% If set to yes, primitive coordinates are used for the polynomials and the right-hand side
302 !% of the Laplacian is computed using the metric tensor and the trace of the Hessian.
303 !% If set to no, Cartesian coordinates are used for the polynomials and the right-hand side
304 !% of the Laplacian reduces to the Cartesian case.
305 !% For some non-orthogonal grids, using primitve coordinates is necessary becasuse the polynomials
306 !% become linearly dependent, thus the corresponding matrix cannot be inverted. For some curvilinear
307 !% coordinate systems, using Cartesian coordinates is more accurate.
308 !%
309 !% By default, use primitve coordinates except for curvilinear meshes.
310 !%End
311 stencil_primitive_coordinates = .not. coord_system%local_basis
312 call parse_variable(namespace, 'StencilPrimitiveCoordinates', stencil_primitive_coordinates, der%stencil_primitive_coordinates)
313
314 ! if needed, der%masses should be initialized in modelmb_particles_init
315 safe_allocate(der%masses(1:space%dim))
316 der%masses = m_one
317
318 ! construct lapl and grad structures
319 safe_allocate(der%op(1:der%dim + 1))
320 der%grad => der%op
321 der%lapl => der%op(der%dim + 1)
322
323 call derivatives_get_stencil_lapl(der, der%lapl, space, coord_system)
325
326 ! find out how many ghost points we need in each dimension
327 safe_allocate(der%n_ghost(1:der%dim))
328 der%n_ghost(:) = 0
329 do idir = 1, der%dim
330 der%n_ghost(idir) = maxval(abs(der%lapl%stencil%points(idir, :)))
331 end do
332
333 nullify(der%coarser)
334 nullify(der%finer)
335 nullify(der%to_coarser)
336 nullify(der%to_finer)
337
338 if (accel_is_enabled()) then
339 ! Check if we can build the uvw_to_xyz kernel
340 select type (coord_system)
341 type is (cartesian_t)
342 ! In this case one does not need to call the kernel, so all is fine
343 class default
344 if (der%dim > 3) then
345 message(1) = "Calculation of derivatives on the GPU with dimension > 3 are only implemented for Cartesian coordinates."
346 call messages_fatal(1, namespace=namespace)
347 else
348 write(flags, '(A,I1.1)') ' -DDIMENSION=', der%dim
349 call accel_kernel_build(kernel_uvw_xyz, 'uvw_to_xyz.cl', 'uvw_to_xyz', flags)
350 end if
351 end select
352 call accel_kernel_build(kernel_dcurl, 'curl.cl', 'dcurl', flags = '-DRTYPE_DOUBLE')
353 call accel_kernel_build(kernel_zcurl, 'curl.cl', 'zcurl', flags = '-DRTYPE_COMPLEX')
354 end if
355
356 pop_sub(derivatives_init)
357 end subroutine derivatives_init
358
359
360 ! ---------------------------------------------------------
361 subroutine derivatives_end(der)
362 type(derivatives_t), intent(inout) :: der
363
364 integer :: idim
365
366 push_sub(derivatives_end)
367
368 assert(allocated(der%op))
369
370 do idim = 1, der%dim+1
371 call nl_operator_end(der%op(idim))
372 end do
373
374 safe_deallocate_a(der%masses)
375
376 safe_deallocate_a(der%n_ghost)
377
378 safe_deallocate_a(der%op)
379 nullify(der%lapl, der%grad)
380
381 nullify(der%coarser)
382 nullify(der%finer)
383 nullify(der%to_coarser)
384 nullify(der%to_finer)
385
386 call boundaries_end(der%boundaries)
387
388 pop_sub(derivatives_end)
389 end subroutine derivatives_end
390
391
392 ! ---------------------------------------------------------
393 subroutine derivatives_get_stencil_lapl(der, lapl, space, coord_system, name, order)
394 type(derivatives_t), intent(in) :: der
395 type(nl_operator_t), intent(inout) :: lapl
396 class(space_t), intent(in) :: space
397 class(coordinate_system_t), intent(in) :: coord_system
398 character(len=80), optional, intent(in) :: name
399 integer, optional, intent(in) :: order
400
401 character(len=80) :: name_
402 integer :: order_
403
405
406 name_ = optional_default(name, "Laplacian")
407 order_ = optional_default(order, der%order)
408
409 ! initialize nl operator
410 call nl_operator_init(lapl, name_)
411
412 ! create stencil
413 select case (der%stencil_type)
414 case (der_star, der_variational)
415 call stencil_star_get_lapl(lapl%stencil, der%dim, order_)
416 case (der_cube)
417 call stencil_cube_get_lapl(lapl%stencil, der%dim, order_)
418 case (der_starplus)
419 call stencil_starplus_get_lapl(lapl%stencil, der%dim, order_)
420 case (der_stargeneral)
421 call stencil_stargeneral_get_arms(lapl%stencil, space%dim, coord_system)
422 call stencil_stargeneral_get_lapl(lapl%stencil, der%dim, order_)
423 end select
424
426 end subroutine derivatives_get_stencil_lapl
427
428
429 ! ---------------------------------------------------------
431 subroutine derivatives_lapl_diag(der, lapl)
432 type(derivatives_t), intent(in) :: der
433 real(real64), intent(out) :: lapl(:)
434
435 push_sub(derivatives_lapl_diag)
436
437 assert(ubound(lapl, dim=1) >= der%mesh%np)
438
439 ! the Laplacian is a real operator
440 call dnl_operator_operate_diag(der%lapl, lapl)
441
442 pop_sub(derivatives_lapl_diag)
443
444 end subroutine derivatives_lapl_diag
445
446
447 ! ---------------------------------------------------------
448 subroutine derivatives_get_stencil_grad(der)
449 type(derivatives_t), intent(inout) :: der
450
451
452 integer :: ii
453 character :: dir_label
456
457 assert(associated(der%grad))
458
459 ! initialize nl operator
460 do ii = 1, der%dim
461 dir_label = ' '
462 if (ii < 5) dir_label = index2axis(ii)
463
464 call nl_operator_init(der%grad(ii), "Gradient "//dir_label)
465
466 ! create stencil
467 select case (der%stencil_type)
468 case (der_star, der_variational)
469 call stencil_star_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order)
470 case (der_cube)
471 call stencil_cube_get_grad(der%grad(ii)%stencil, der%dim, der%order)
472 case (der_starplus)
473 call stencil_starplus_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order)
474 case (der_stargeneral)
475 ! use the simple star stencil
476 call stencil_star_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order)
477 end select
478 end do
479
481
482 end subroutine derivatives_get_stencil_grad
483
484 ! ---------------------------------------------------------
530 subroutine derivatives_build(der, namespace, space, mesh, qvector, regenerate, verbose)
531 type(derivatives_t), intent(inout) :: der
532 type(namespace_t), intent(in) :: namespace
533 class(space_t), intent(in) :: space
534 class(mesh_t), target, intent(in) :: mesh
535 real(real64), optional, intent(in) :: qvector(:)
537 logical, optional, intent(in) :: regenerate
538 logical, optional, intent(in) :: verbose
539
540 integer :: i
541 logical :: const_w_
542 integer :: np_zero_bc
543
544 push_sub(derivatives_build)
545
546 if (.not. optional_default(regenerate, .false.)) then
547 call boundaries_init(der%boundaries, namespace, space, mesh, qvector)
548 end if
549
550 assert(allocated(der%op))
551 assert(der%stencil_type >= der_star .and. der%stencil_type <= der_stargeneral)
552 assert(.not.(der%stencil_type == der_variational .and. mesh%use_curvilinear))
553
554 der%mesh => mesh ! make a pointer to the underlying mesh
555
556 const_w_ = .true.
557
558 ! need non-constant weights for curvilinear and scattering meshes
559 if (mesh%use_curvilinear) const_w_ = .false.
560
561 np_zero_bc = 0
562
563 ! build operators
564 do i = 1, der%dim+1
565 if (optional_default(regenerate, .false.)) then
566 safe_deallocate_a(der%op(i)%w)
567 end if
568 call nl_operator_build(space, mesh, der%op(i), der%mesh%np, const_w = const_w_, &
569 regenerate=regenerate)
570 np_zero_bc = max(np_zero_bc, nl_operator_np_zero_bc(der%op(i)))
571 end do
572
573 assert(np_zero_bc > mesh%np .and. np_zero_bc <= mesh%np_part)
574
575 select case (der%stencil_type)
576
577 case (der_star, der_starplus, der_stargeneral) ! Laplacian and gradient have different stencils
578 do i = 1, der%dim + 1
579
580 call derivatives_make_discretization(namespace, der, 1, i, der%op(i:i), verbose=verbose)
581 end do
582
583 case (der_cube)
584 ! gradients have the same stencils, so use one call to derivatives_make_discretization
585 ! to solve the linear equation once for several right-hand sides
586 call derivatives_make_discretization(namespace, der, der%dim, der%dim, der%op(:), verbose=verbose)
587 ! the polynomials are evaluated differently for the Laplacian
588 call derivatives_make_discretization(namespace, der, 1, der%dim+1, der%op(der%dim+1:der%dim+1), verbose=verbose)
589
590 case (der_variational)
591 ! we have the explicit coefficients
592 call stencil_variational_coeff_lapl(der%dim, der%order, mesh%spacing, der%lapl, alpha = der%lapl_cutoff)
593 end select
594
595 pop_sub(derivatives_build)
596
597 end subroutine derivatives_build
598
599 ! ---------------------------------------------------------
600 subroutine stencil_pol_grad(stencil_type, dim, direction, order, polynomials)
601 integer, intent(in) :: stencil_type
602 integer, intent(in) :: dim
603 integer, intent(in) :: direction
604 integer, intent(in) :: order
605 integer, intent(inout) :: polynomials(:, :)
606
607 select case (stencil_type)
608 case (der_star, der_stargeneral)
609 call stencil_star_polynomials_grad(direction, order, polynomials)
610 case (der_starplus)
611 call stencil_starplus_pol_grad(dim, direction, order, polynomials)
612 case (der_cube)
613 ! same stencil for gradient and Laplacian in this case
614 call stencil_cube_polynomials_lapl(dim, order, polynomials)
615 end select
616 end subroutine stencil_pol_grad
617
618 ! ---------------------------------------------------------
619 subroutine stencil_pol_lapl(stencil_type, stencil, dim, order, polynomials)
620 integer, intent(in) :: stencil_type
621 type(stencil_t), intent(in) :: stencil
622 integer, intent(in) :: dim
623 integer, intent(in) :: order
624 integer, intent(inout) :: polynomials(:, :)
625
626 select case (stencil_type)
627 case (der_star)
628 call stencil_star_polynomials_lapl(dim, order, polynomials)
629 case (der_starplus)
630 call stencil_starplus_pol_lapl(dim, order, polynomials)
631 case (der_stargeneral)
632 call stencil_stargeneral_pol_lapl(stencil, dim, order, polynomials)
633 case (der_cube)
634 call stencil_cube_polynomials_lapl(dim, order, polynomials)
635 end select
636 end subroutine stencil_pol_lapl
637
638 ! ---------------------------------------------------------
639 subroutine get_rhs_lapl(coord_system, polynomials, x, rhs, stencil_primitive_coordinates)
640 class(coordinate_system_t), intent(in) :: coord_system
641 integer, intent(in) :: polynomials(:,:)
642 real(real64), intent(in) :: x(:)
643 real(real64), intent(out) :: rhs(:)
644 logical, intent(in) :: stencil_primitive_coordinates
645
646 integer :: i, j, k, dim
647 real(real64) :: metric(1:size(polynomials, dim=1), 1:size(polynomials, dim=1)), chi(1:size(polynomials, dim=1))
648 real(real64) :: jacobian(1:size(polynomials, dim=1), 1:size(polynomials, dim=1))
649 real(real64) :: hessian_trace(1:size(polynomials, dim=1))
650 integer :: powers(0:2)
651
652 push_sub(get_rhs_lapl)
653
654 dim = size(polynomials, dim=1)
655
656 if (stencil_primitive_coordinates) then
657 ! metric and trace of Hessian in primitive coordinates
658 select type (coord => coord_system)
659 class is (affine_coordinates_t)
660 metric(1:dim, 1:dim) = matmul(coord%basis%change_of_basis_matrix(1:dim, 1:dim), &
661 transpose(coord%basis%change_of_basis_matrix(1:dim, 1:dim)))
662 hessian_trace = m_zero
663 class is (curv_gygi_t)
664 call curv_gygi_jacobian(coord, x, chi, jacobian)
665 metric(1:dim, 1:dim) = matmul(jacobian(1:dim, 1:dim), transpose(jacobian(1:dim, 1:dim)))
666 call curv_gygi_hessian_trace(coord, x, hessian_trace)
667 class default
668 message(1) = "Weight computation not implemented for the coordinate system chosen."
669 call messages_fatal(1)
670 end select
671 else
672 ! metric and trace of Hessian in Cartesian coordinates
673 metric = m_zero
674 do i = 1, dim
675 metric(i, i) = m_one
676 end do
677 hessian_trace = m_zero
678 end if
679
680 ! find right-hand side for operator
681 rhs(:) = m_zero
682 do j = 1, size(polynomials, dim=2)
683 ! count the powers of the polynomials
684 powers = 0
685 do i = 1, dim
686 if (polynomials(i, j) <= 2) then
687 powers(polynomials(i, j)) = powers(polynomials(i, j)) + 1
688 end if
689 end do
690
691 ! find all polynomials for which exactly one term is quadratic
692 ! for these, the Laplacian on the polynomial is 2*metric(i, i)
693 if (powers(2) == 1 .and. powers(0) == dim - 1) then
694 do i = 1, dim
695 if (polynomials(i, j) == 2) then
696 rhs(j) = m_two*metric(i, i)
697 end if
698 end do
699 end if
700 ! find all polynomials for which exactly two terms are linear
701 ! for these, the Laplacian on the polynomial is metric(i, k) + metric(k, i)
702 if (powers(1) == 2 .and. powers(0) == dim - 2) then
703 do i = 1, dim
704 if (polynomials(i, j) == 1) then
705 do k = i+1, dim
706 if (polynomials(k, j) == 1) then
707 rhs(j) = metric(i, k) + metric(k, i)
708 end if
709 end do
710 end if
711 end do
712 end if
713 ! find all polynomials for which exactly one term is linear
714 ! for these, the Laplacian on the polynomial is hessian_trace(i)
715 ! this term is only non-zero for curvilinear coordinates with a varying metric tensor
716 if (powers(1) == 1 .and. powers(0) == dim - 1) then
717 do i = 1, dim
718 if (polynomials(i, j) == 1) then
719 rhs(j) = hessian_trace(i)
720 end if
721 end do
722 end if
723 end do
724
725 pop_sub(get_rhs_lapl)
726 end subroutine get_rhs_lapl
727
728 ! ---------------------------------------------------------
729 subroutine get_rhs_grad(polynomials, dir, rhs)
730 integer, intent(in) :: polynomials(:,:)
731 integer, intent(in) :: dir
732 real(real64), intent(out) :: rhs(:)
733
734 integer :: j, k, dim
735 logical :: this_one
736
737 push_sub(get_rhs_grad)
738
739 dim = size(polynomials, dim=1)
740
741 ! find right-hand side for operator
742 rhs(:) = m_zero
743 do j = 1, size(polynomials, dim=2)
744 this_one = .true.
745 do k = 1, dim
746 if (k == dir .and. polynomials(k, j) /= 1) this_one = .false.
747 if (k /= dir .and. polynomials(k, j) /= 0) this_one = .false.
748 end do
749 if (this_one) rhs(j) = m_one
750 end do
751
752 pop_sub(get_rhs_grad)
753 end subroutine get_rhs_grad
754
755
756 ! ---------------------------------------------------------
757 subroutine derivatives_make_discretization(namespace, der, nderiv, ideriv, op, name, verbose, order)
758 type(namespace_t), intent(in) :: namespace
759 type(derivatives_t), intent(in) :: der
760 integer, intent(in) :: nderiv
761 integer, intent(in) :: ideriv
762 type(nl_operator_t), intent(inout) :: op(:)
763 character(len=80), optional, intent(in) :: name
764 logical, optional, intent(in) :: verbose
765 integer, optional, intent(in) :: order
766
767 integer :: p, p_max, i, j, k, pow_max, order_
768 real(real64) :: x(der%dim)
769 real(real64), allocatable :: mat(:,:), sol(:,:), powers(:,:)
770 integer, allocatable :: pol(:,:)
771 real(real64), allocatable :: rhs(:,:)
772 character(len=80) :: name_
773
775
776 order_ = optional_default(order, der%order)
777
778 safe_allocate(rhs(1:op(1)%stencil%size, nderiv))
779 safe_allocate(pol(1:der%dim, 1:op(1)%stencil%size))
780 ! get polynomials
781 if (nderiv == 1) then
782 if (ideriv <= der%dim) then ! gradient
783 call stencil_pol_grad(der%stencil_type, der%dim, ideriv, order_, pol)
784 name_ = index2axis(ideriv) // "-gradient"
785 else ! Laplacian
786 call stencil_pol_lapl(der%stencil_type, op(1)%stencil, der%dim, order_, pol)
787 name_ = "Laplacian"
788 end if
789 else if (nderiv == der%dim) then
790 ! in this case, we compute the weights for the gradients at once
791 ! because they have the same stencils
792 ! this is only used for the cube stencil
793 name_ = "gradients"
794 call stencil_pol_grad(der%stencil_type, der%dim, 1, order_, pol)
795 else
796 message(1) = "Error: derivatives_make_discretization can only be called with nderiv = 1 or der%dim."
797 call messages_fatal(1)
798 end if
799
800 safe_allocate(mat(1:op(1)%stencil%size, 1:op(1)%stencil%size))
801 safe_allocate(sol(1:op(1)%stencil%size, 1:nderiv))
802
803 if (optional_default(verbose, .true.)) then
804 name_ = optional_default(name, name_)
805 message(1) = 'Info: Generating weights for finite-difference discretization of ' // trim(name_)
806 call messages_info(1, namespace=namespace)
807 end if
808
809 ! use to generate power lookup table
810 pow_max = maxval(pol)
811 safe_allocate(powers(1:der%dim, 0:pow_max))
812 powers(:,:) = m_zero
813 powers(:,0) = m_one
814
815 p_max = op(1)%np
816 if (op(1)%const_w) p_max = 1
817
818 do p = 1, p_max
819 ! get polynomials and right-hand side
820 ! we need the current position for the right-hand side for curvilinear coordinates
821 x = der%mesh%x(p, :)
822 if (nderiv == 1) then
823 if (ideriv <= der%dim) then ! gradient
824 call get_rhs_grad(pol, ideriv, rhs(:,1))
825 name_ = index2axis(ideriv) // "-gradient"
826 else ! Laplacian
827 call get_rhs_lapl(der%mesh%coord_system, pol, x, rhs(:,1), der%stencil_primitive_coordinates)
828 end if
829 else if (nderiv == der%dim) then
830 ! in this case, we compute the weights for the gradients at once
831 ! because they have the same stencils
832 do i = 1, der%dim
833 call get_rhs_grad(pol, i, rhs(:,i))
834 end do
835 end if
836
837 ! first polynomial is just a constant
838 mat(1,:) = m_one
839 ! i indexes the point in the stencil
840 do i = 1, op(1)%stencil%size
841 if (ideriv <= der%dim) then
842 ! for gradients, we always need to evaluate the polynomials in primitive coordinates
843 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
844 else
845 ! for the Laplacian, we can evaluate the polynomials in primitive or Cartesian coordinates
846 if (der%stencil_primitive_coordinates) then
847 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
848 else
849 x = der%mesh%x(p + op(1)%ri(i, op(1)%rimap(p)), :) - der%mesh%x(p, :)
850 end if
851 end if
852
853 ! NB: these masses are applied on the cartesian directions. Should add a check for non-orthogonal axes
854 x = x*sqrt(der%masses)
855
856 ! calculate powers
857 powers(:, 1) = x
858 do k = 2, pow_max
859 powers(:, k) = x*powers(:, k-1)
860 end do
861
862 ! generate the matrix
863 ! j indexes the polynomial being used
864 do j = 2, op(1)%stencil%size
865 mat(j, i) = powers(1, pol(1, j))
866 do k = 2, der%dim
867 mat(j, i) = mat(j, i)*powers(k, pol(k, j))
868 end do
869 end do
870 end do ! loop over i = point in stencil
871
872 ! linear problem to solve for derivative weights:
873 ! mat * sol = rhs
874 call lalg_linsyssolve(op(1)%stencil%size, nderiv, mat, rhs, sol)
875
876 ! for the cube stencil, all derivatives are calculated at once, so assign
877 ! the correct solution to each operator
878 do i = 1, nderiv
879 op(i)%w(:, p) = sol(:, i)
880 end do
881
882 end do ! loop over points p
883
884 do i = 1, nderiv
886 end do
887
888 ! In case of constant weights, we store the weights of the Laplacian on the GPU, as this
889 ! saves many unecessary transfers
890 do i = 1, nderiv
893 end do
894
895
896 safe_deallocate_a(mat)
897 safe_deallocate_a(sol)
898 safe_deallocate_a(powers)
899 safe_deallocate_a(pol)
900 safe_deallocate_a(rhs)
901
904
905#ifdef HAVE_MPI
906 ! ---------------------------------------------------------
907 logical function derivatives_overlap(this) result(overlap)
908 type(derivatives_t), intent(in) :: this
909
910 push_sub(derivatives_overlap)
911
912 overlap = this%comm_method /= blocking
913
914 pop_sub(derivatives_overlap)
915 end function derivatives_overlap
916#endif
917
918 ! ---------------------------------------------------------
919 subroutine derivatives_get_lapl(this, namespace, op, space, name, order)
920 type(derivatives_t), intent(in) :: this
921 type(namespace_t), intent(in) :: namespace
922 type(nl_operator_t), intent(inout) :: op(:)
923 class(space_t), intent(in) :: space
924 character(len=80), intent(in) :: name
925 integer, intent(in) :: order
926
927 push_sub(derivatives_get_lapl)
928
929 call derivatives_get_stencil_lapl(this, op(1), space, this%mesh%coord_system, name, order)
930 call nl_operator_build(space, this%mesh, op(1), this%mesh%np, const_w = .not. this%mesh%use_curvilinear)
931 call derivatives_make_discretization(namespace, this, 1, this%dim+1, op(1:1), name, order=order)
932
933 pop_sub(derivatives_get_lapl)
934 end subroutine derivatives_get_lapl
935
936 ! ---------------------------------------------------------
949 function derivatives_get_inner_boundary_mask(this) result(mask)
950 type(derivatives_t), intent(in) :: this
951
952 logical :: mask(1:this%mesh%np)
953 integer :: ip, is, index
954
955 mask = .false.
956 ! Loop through all points in the grid
957 do ip = 1, this%mesh%np
958 ! For each of them, loop through all points in the stencil
959 do is = 1, this%lapl%stencil%size
960 ! Get the index of the point obtained as: grid_point + displament_due_to_stencil
961 index = nl_operator_get_index(this%lapl, is, ip)
962 ! Check whether the displaced point if outsude the grid. Is so, it belongs to the mask
963 if (index > this%mesh%np + this%mesh%pv%np_ghost) then
964 mask(ip) = .true.
965 exit
966 end if
967 end do
968 end do
969
971
972 ! ---------------------------------------------------------
977 real(real64) function derivatives_lapl_get_max_eigenvalue(this)
978 type(derivatives_t), intent(in) :: this
979
980 integer :: i
981
982 push_sub(derivatives_lapl_get_max_eigenvalue)
983
984 derivatives_lapl_get_max_eigenvalue = m_zero
985 if ((this%stencil_type == der_star .or. this%stencil_type == der_stargeneral) &
986 .and. .not. this%mesh%use_curvilinear) then
987 ! use Fourier transform of stencil evaluated at the maximum phase
988 do i = 1, this%lapl%stencil%size
989 derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
990 (-1)**maxval(abs(this%lapl%stencil%points(:, i)))*this%lapl%w(i, 1)
991 end do
992 derivatives_lapl_get_max_eigenvalue = abs(derivatives_lapl_get_max_eigenvalue)
993 else
994 ! use upper bound from continuum for other stencils
995 do i = 1, this%dim
996 derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
997 m_pi**2/this%mesh%spacing(i)**2
998 end do
999 end if
1001 pop_sub(derivatives_lapl_get_max_eigenvalue)
1003
1004
1005#include "undef.F90"
1006#include "real.F90"
1007#include "derivatives_inc.F90"
1008
1009#include "undef.F90"
1010#include "complex.F90"
1011#include "derivatives_inc.F90"
1013end module derivatives_oct_m
1014
1015!! Local Variables:
1016!! mode: f90
1017!! coding: utf-8
1018!! End:
double sqrt(double __x) __attribute__((__nothrow__
subroutine, public accel_kernel_build(this, file_name, kernel_name, flags)
Definition: accel.F90:2067
pure logical function, public accel_is_enabled()
Definition: accel.F90:402
This module implements batches of mesh functions.
Definition: batch.F90:133
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:116
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:122
subroutine, public boundaries_end(this)
Definition: boundaries.F90:428
subroutine, public boundaries_init(this, namespace, space, mesh, qvector)
initialize the boundary contitions
Definition: boundaries.F90:227
This module implements the curvilinear coordinates given in E.L. Briggs, D.J. Sullivan,...
This module implements the curvilinear coordinates given in F. Gygi and G. Galli, PRB 52 R2229 (1996)...
Definition: curv_gygi.F90:117
pure subroutine, public curv_gygi_jacobian(this, xx, chi, jac, natoms)
Definition: curv_gygi.F90:445
pure subroutine, public curv_gygi_hessian_trace(this, xx, hessian_trace, natoms)
Definition: curv_gygi.F90:527
This module implements the curvilinear coordinates given in N. A. Modine, G. Zumbach,...
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_test(this, namespace, repetitions, min_blocksize, max_blocksize)
unit test for derivatives
subroutine, public dderivatives_batch_grad(der, ffb, opffb, ghost_update, set_bc, to_cartesian, metric, factor)
apply the gradient to a batch of mesh functions
subroutine, public zderivatives_curl(der, ff, op_ff, ghost_update, set_bc)
apply the curl operator to a vector of mesh functions
subroutine, public dderivatives_perform(op, der, ff, op_ff, ghost_update, set_bc, factor)
apply a nl_operator to a mesh function
subroutine, public derivatives_lapl_diag(der, lapl)
Returns the diagonal elements of the Laplacian, needed for preconditioning.
integer, parameter, public der_cube
subroutine, public dderivatives_batch_curl_from_gradient(der, ffb, gradb)
calculate the curl from a batch and its gradient
subroutine, public derivatives_build(der, namespace, space, mesh, qvector, regenerate, verbose)
build the derivatives object:
subroutine stencil_pol_grad(stencil_type, dim, direction, order, polynomials)
subroutine stencil_pol_lapl(stencil_type, stencil, dim, order, polynomials)
subroutine, public zderivatives_batch_curl_from_gradient(der, ffb, gradb)
calculate the curl from a batch and its gradient
subroutine, public zderivatives_batch_div(der, ffb, opffb, ghost_update, set_bc, to_cartesian)
calculate the divergence of a vector of batches
subroutine, public dderivatives_partial(der, ff, op_ff, dir, ghost_update, set_bc)
apply the partial derivative along dir to a mesh function
subroutine, public dderivatives_batch_perform(op, der, ff, opff, ghost_update, set_bc, factor, async)
apply an operator to a bach of mesh functions
subroutine, public zderivatives_div(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the divergence operator to a vector of mesh functions
integer, parameter der_bc_period
boundary is periodic
subroutine, public dderivatives_batch_finish(handle, async)
apply a non-local operator to a batch (2nd part)
subroutine, public dderivatives_curl(der, ff, op_ff, ghost_update, set_bc)
apply the curl operator to a vector of mesh functions
subroutine derivatives_get_stencil_grad(der)
subroutine, public dderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine derivatives_get_stencil_lapl(der, lapl, space, coord_system, name, order)
subroutine, public derivatives_init(der, namespace, space, coord_system, order)
subroutine, public derivatives_end(der)
integer, parameter der_bc_zero_df
first derivative of the function is zero
subroutine get_rhs_lapl(coord_system, polynomials, x, rhs, stencil_primitive_coordinates)
subroutine, public dderivatives_batch_div(der, ffb, opffb, ghost_update, set_bc, to_cartesian)
calculate the divergence of a vector of batches
subroutine, public derivatives_get_lapl(this, namespace, op, space, name, order)
subroutine, public zderivatives_batch_perform(op, der, ff, opff, ghost_update, set_bc, factor, async)
apply an operator to a bach of mesh functions
subroutine, public zderivatives_batch_grad(der, ffb, opffb, ghost_update, set_bc, to_cartesian, metric, factor)
apply the gradient to a batch of mesh functions
subroutine, public dderivatives_batch_start(op, der, ff, opff, handle, ghost_update, set_bc, factor)
apply a non-local operator to a batch (1st part)
subroutine get_rhs_grad(polynomials, dir, rhs)
subroutine, public dderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
real(real64) function, public derivatives_lapl_get_max_eigenvalue(this)
Get maximum eigenvalue of discrete Laplacian. For the star and star_general stencils,...
logical function derivatives_overlap(this)
subroutine, public dderivatives_batch_curl(der, ffb, ghost_update, set_bc)
apply the curl to a batch of mesh functions
subroutine, public zderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine, public zderivatives_partial(der, ff, op_ff, dir, ghost_update, set_bc)
apply the partial derivative along dir to a mesh function
subroutine, public zderivatives_batch_start(op, der, ff, opff, handle, ghost_update, set_bc, factor)
apply a non-local operator to a batch (1st part)
logical function, dimension(1:this%mesh%np), public derivatives_get_inner_boundary_mask(this)
This function tells whether a point in the grid is contained in a layer of the width of the stencil b...
integer, parameter, public der_starplus
subroutine derivatives_make_discretization(namespace, der, nderiv, ideriv, op, name, verbose, order)
subroutine, public zderivatives_perform(op, der, ff, op_ff, ghost_update, set_bc, factor)
apply a nl_operator to a mesh function
subroutine, public zderivatives_batch_curl(der, ffb, ghost_update, set_bc)
apply the curl to a batch of mesh functions
subroutine, public zderivatives_test(this, namespace, repetitions, min_blocksize, max_blocksize)
unit test for derivatives
integer, parameter, public der_variational
subroutine, public zderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
subroutine, public zderivatives_batch_finish(handle, async)
apply a non-local operator to a batch (2nd part)
subroutine, public dderivatives_div(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the divergence operator to a vector of mesh functions
integer, parameter non_blocking
integer, parameter, public der_stargeneral
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
real(real64), parameter, public m_one
Definition: global.F90:189
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
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_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
This module defines non-local operators.
subroutine, public dnl_operator_operate_diag(op, fo)
subroutine, public nl_operator_init(op, label)
initialize an instance of a non-local operator by setting the label
subroutine, public nl_operator_update_gpu_buffers(op)
subroutine, public nl_operator_output_weights(this)
subroutine, public nl_operator_end(op)
subroutine, public nl_operator_build(space, mesh, op, np, const_w, regenerate)
integer pure function, public nl_operator_np_zero_bc(op)
integer pure function, public nl_operator_get_index(op, is, ip)
subroutine, public nl_operator_allocate_gpu_buffers(op)
Some general things and nomenclature:
Definition: par_vec.F90:171
This module defines routines, generating operators for a cubic stencil.
subroutine, public stencil_cube_get_lapl(this, dim, order)
subroutine, public stencil_cube_polynomials_lapl(dim, order, pol)
subroutine, public stencil_cube_get_grad(this, dim, order)
This module defines stencils used in Octopus.
Definition: stencil.F90:135
This module defines routines, generating operators for a star stencil.
subroutine, public stencil_star_get_grad(this, dim, dir, order)
subroutine, public stencil_star_get_lapl(this, dim, order)
subroutine, public stencil_star_polynomials_grad(dir, order, pol)
subroutine, public stencil_star_polynomials_lapl(dim, order, pol)
This module defines routines, generating operators for a generalized star stencil.
subroutine, public stencil_stargeneral_get_arms(this, dim, coord_system)
Finds the direction of the arms of the star-general stencil as described in Natan et al....
subroutine, public stencil_stargeneral_pol_lapl(this, dim, order, pol)
subroutine, public stencil_stargeneral_get_lapl(this, dim, order)
This module defines routines, generating operators for a stencil consisting of a star and a cross....
subroutine, public stencil_starplus_pol_grad(dim, dir, order, pol)
subroutine, public stencil_starplus_pol_lapl(dim, order, pol)
subroutine, public stencil_starplus_get_lapl(this, dim, order)
subroutine, public stencil_starplus_get_grad(this, dim, dir, order)
Implements the variational discretization of the Laplacian as proposed by P. Maragakis,...
subroutine, public stencil_variational_coeff_lapl(dim, order, h, lapl, alpha)
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
character pure function, public index2axis(idir)
Definition: utils.F90:202
abstract class to describe coordinate systems
handle to transfer data from the start() to finish() calls.
class representing derivatives
Describes mesh distribution to nodes.
Definition: mesh.F90:186
data type for non local operators
The class representing the stencil, which is used for non-local mesh operations.
Definition: stencil.F90:163
int true(void)