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, chi, rhs, stencil_primitive_coordinates)
640 class(coordinate_system_t), intent(in) :: coord_system
641 integer, intent(in) :: polynomials(:,:)
642 real(real64), intent(in) :: chi(:)
643 real(real64), intent(out) :: rhs(:)
644 logical, intent(in) :: stencil_primitive_coordinates
645
646 integer :: i, j, k, dim
647 real(real64) :: metric_inverse(1:size(polynomials, dim=1), 1:size(polynomials, dim=1))
648 real(real64) :: hessian_trace(1:size(polynomials, dim=1))
649 integer :: powers(0:2)
650
651 push_sub(get_rhs_lapl)
652
653 dim = size(polynomials, dim=1)
654
655 if (stencil_primitive_coordinates) then
656 ! inverse metric and trace of Hessian in primitive coordinates
657 metric_inverse = coord_system%metric_inverse(chi)
658 hessian_trace = coord_system%trace_hessian(chi)
659 else
660 ! inverse metric and trace of Hessian in Cartesian coordinates
661 metric_inverse = m_zero
662 do i = 1, dim
663 metric_inverse(i, i) = m_one
664 end do
665 hessian_trace = m_zero
666 end if
667
668 ! find right-hand side for operator
669 rhs(:) = m_zero
670 do j = 1, size(polynomials, dim=2)
671 ! count the powers of the polynomials
672 powers = 0
673 do i = 1, dim
674 if (polynomials(i, j) <= 2) then
675 powers(polynomials(i, j)) = powers(polynomials(i, j)) + 1
676 end if
677 end do
678
679 ! find all polynomials for which exactly one term is quadratic
680 ! for these, the Laplacian on the polynomial is 2*metric_inverse(i, i)
681 if (powers(2) == 1 .and. powers(0) == dim - 1) then
682 do i = 1, dim
683 if (polynomials(i, j) == 2) then
684 rhs(j) = m_two*metric_inverse(i, i)
685 end if
686 end do
687 end if
688 ! find all polynomials for which exactly two terms are linear
689 ! for these, the Laplacian on the polynomial is metric_inverse(i, k) + metric_inverse(k, i)
690 if (powers(1) == 2 .and. powers(0) == dim - 2) then
691 do i = 1, dim
692 if (polynomials(i, j) == 1) then
693 do k = i+1, dim
694 if (polynomials(k, j) == 1) then
695 rhs(j) = metric_inverse(i, k) + metric_inverse(k, i)
696 end if
697 end do
698 end if
699 end do
700 end if
701 ! find all polynomials for which exactly one term is linear
702 ! for these, the Laplacian on the polynomial is hessian_trace(i)
703 ! this term is only non-zero for curvilinear coordinates with a varying metric tensor
704 if (powers(1) == 1 .and. powers(0) == dim - 1) then
705 do i = 1, dim
706 if (polynomials(i, j) == 1) then
707 rhs(j) = hessian_trace(i)
708 end if
709 end do
710 end if
711 end do
713 pop_sub(get_rhs_lapl)
714 end subroutine get_rhs_lapl
715
716 ! ---------------------------------------------------------
717 subroutine get_rhs_grad(polynomials, dir, rhs)
718 integer, intent(in) :: polynomials(:,:)
719 integer, intent(in) :: dir
720 real(real64), intent(out) :: rhs(:)
721
722 integer :: j, k, dim
723 logical :: this_one
724
725 push_sub(get_rhs_grad)
726
727 dim = size(polynomials, dim=1)
728
729 ! find right-hand side for operator
730 rhs(:) = m_zero
731 do j = 1, size(polynomials, dim=2)
732 this_one = .true.
733 do k = 1, dim
734 if (k == dir .and. polynomials(k, j) /= 1) this_one = .false.
735 if (k /= dir .and. polynomials(k, j) /= 0) this_one = .false.
736 end do
737 if (this_one) rhs(j) = m_one
738 end do
739
740 pop_sub(get_rhs_grad)
741 end subroutine get_rhs_grad
742
743
744 ! ---------------------------------------------------------
745 subroutine derivatives_make_discretization(namespace, der, nderiv, ideriv, op, name, verbose, order)
746 type(namespace_t), intent(in) :: namespace
747 type(derivatives_t), intent(in) :: der
748 integer, intent(in) :: nderiv
749 integer, intent(in) :: ideriv
750 type(nl_operator_t), intent(inout) :: op(:)
751 character(len=80), optional, intent(in) :: name
752 logical, optional, intent(in) :: verbose
753 integer, optional, intent(in) :: order
754
755 integer :: p, p_max, i, j, k, pow_max, order_
756 real(real64) :: x(der%dim)
757 real(real64), allocatable :: mat(:,:), sol(:,:), powers(:,:)
758 integer, allocatable :: pol(:,:)
759 real(real64), allocatable :: rhs(:,:)
760 character(len=80) :: name_
761
763
764 order_ = optional_default(order, der%order)
765
766 safe_allocate(rhs(1:op(1)%stencil%size, nderiv))
767 safe_allocate(pol(1:der%dim, 1:op(1)%stencil%size))
768 ! get polynomials
769 if (nderiv == 1) then
770 if (ideriv <= der%dim) then ! gradient
771 call stencil_pol_grad(der%stencil_type, der%dim, ideriv, order_, pol)
772 name_ = index2axis(ideriv) // "-gradient"
773 else ! Laplacian
774 call stencil_pol_lapl(der%stencil_type, op(1)%stencil, der%dim, order_, pol)
775 name_ = "Laplacian"
776 end if
777 else if (nderiv == der%dim) then
778 ! in this case, we compute the weights for the gradients at once
779 ! because they have the same stencils
780 ! this is only used for the cube stencil
781 name_ = "gradients"
782 call stencil_pol_grad(der%stencil_type, der%dim, 1, order_, pol)
783 else
784 message(1) = "Error: derivatives_make_discretization can only be called with nderiv = 1 or der%dim."
785 call messages_fatal(1)
786 end if
787
788 safe_allocate(mat(1:op(1)%stencil%size, 1:op(1)%stencil%size))
789 safe_allocate(sol(1:op(1)%stencil%size, 1:nderiv))
790
791 if (optional_default(verbose, .true.)) then
792 name_ = optional_default(name, name_)
793 message(1) = 'Info: Generating weights for finite-difference discretization of ' // trim(name_)
794 call messages_info(1, namespace=namespace)
795 end if
796
797 ! use to generate power lookup table
798 pow_max = maxval(pol)
799 safe_allocate(powers(1:der%dim, 0:pow_max))
800 powers(:,:) = m_zero
801 powers(:,0) = m_one
802
803 p_max = op(1)%np
804 if (op(1)%const_w) p_max = 1
805
806 do p = 1, p_max
807 ! get polynomials and right-hand side
808 ! we need the current position for the right-hand side for curvilinear coordinates
809 if (nderiv == 1) then
810 if (ideriv <= der%dim) then ! gradient
811 call get_rhs_grad(pol, ideriv, rhs(:,1))
812 name_ = index2axis(ideriv) // "-gradient"
813 else ! Laplacian
814 call get_rhs_lapl(der%mesh%coord_system, pol, der%mesh%chi(:, p), rhs(:,1), der%stencil_primitive_coordinates)
815 end if
816 else if (nderiv == der%dim) then
817 ! in this case, we compute the weights for the gradients at once
818 ! because they have the same stencils
819 do i = 1, der%dim
820 call get_rhs_grad(pol, i, rhs(:,i))
821 end do
822 end if
823
824 ! first polynomial is just a constant
825 mat(1,:) = m_one
826 ! i indexes the point in the stencil
827 do i = 1, op(1)%stencil%size
828 if (ideriv <= der%dim) then
829 ! for gradients, we always need to evaluate the polynomials in primitive coordinates
830 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
831 else
832 ! for the Laplacian, we can evaluate the polynomials in primitive or Cartesian coordinates
833 if (der%stencil_primitive_coordinates) then
834 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
835 else
836 x = der%mesh%x(p + op(1)%ri(i, op(1)%rimap(p)), :) - der%mesh%x(p, :)
837 end if
838 end if
839
840 ! NB: these masses are applied on the cartesian directions. Should add a check for non-orthogonal axes
841 x = x*sqrt(der%masses)
842
843 ! calculate powers
844 powers(:, 1) = x
845 do k = 2, pow_max
846 powers(:, k) = x*powers(:, k-1)
847 end do
848
849 ! generate the matrix
850 ! j indexes the polynomial being used
851 do j = 2, op(1)%stencil%size
852 mat(j, i) = powers(1, pol(1, j))
853 do k = 2, der%dim
854 mat(j, i) = mat(j, i)*powers(k, pol(k, j))
855 end do
856 end do
857 end do ! loop over i = point in stencil
858
859 ! linear problem to solve for derivative weights:
860 ! mat * sol = rhs
861 call lalg_linsyssolve(op(1)%stencil%size, nderiv, mat, rhs, sol)
862
863 ! for the cube stencil, all derivatives are calculated at once, so assign
864 ! the correct solution to each operator
865 do i = 1, nderiv
866 op(i)%w(:, p) = sol(:, i)
867 end do
868
869 end do ! loop over points p
870
871 do i = 1, nderiv
873 end do
874
875 ! In case of constant weights, we store the weights of the Laplacian on the GPU, as this
876 ! saves many unecessary transfers
877 do i = 1, nderiv
880 end do
881
882
883 safe_deallocate_a(mat)
884 safe_deallocate_a(sol)
885 safe_deallocate_a(powers)
886 safe_deallocate_a(pol)
887 safe_deallocate_a(rhs)
888
891
892#ifdef HAVE_MPI
893 ! ---------------------------------------------------------
894 logical function derivatives_overlap(this) result(overlap)
895 type(derivatives_t), intent(in) :: this
896
897 push_sub(derivatives_overlap)
898
899 overlap = this%comm_method /= blocking
900
901 pop_sub(derivatives_overlap)
902 end function derivatives_overlap
903#endif
904
905 ! ---------------------------------------------------------
906 subroutine derivatives_get_lapl(this, namespace, op, space, name, order)
907 type(derivatives_t), intent(in) :: this
908 type(namespace_t), intent(in) :: namespace
909 type(nl_operator_t), intent(inout) :: op(:)
910 class(space_t), intent(in) :: space
911 character(len=80), intent(in) :: name
912 integer, intent(in) :: order
913
914 push_sub(derivatives_get_lapl)
915
916 call derivatives_get_stencil_lapl(this, op(1), space, this%mesh%coord_system, name, order)
917 call nl_operator_build(space, this%mesh, op(1), this%mesh%np, const_w = .not. this%mesh%use_curvilinear)
918 call derivatives_make_discretization(namespace, this, 1, this%dim+1, op(1:1), name, order=order)
919
920 pop_sub(derivatives_get_lapl)
921 end subroutine derivatives_get_lapl
922
923 ! ---------------------------------------------------------
936 function derivatives_get_inner_boundary_mask(this) result(mask)
937 type(derivatives_t), intent(in) :: this
938
939 logical :: mask(1:this%mesh%np)
940 integer :: ip, is, index
941
942 mask = .false.
943 ! Loop through all points in the grid
944 do ip = 1, this%mesh%np
945 ! For each of them, loop through all points in the stencil
946 do is = 1, this%lapl%stencil%size
947 ! Get the index of the point obtained as: grid_point + displament_due_to_stencil
948 index = nl_operator_get_index(this%lapl, is, ip)
949 ! Check whether the displaced point if outsude the grid. Is so, it belongs to the mask
950 if (index > this%mesh%np + this%mesh%pv%np_ghost) then
951 mask(ip) = .true.
952 exit
953 end if
954 end do
955 end do
956
958
959 ! ---------------------------------------------------------
964 real(real64) function derivatives_lapl_get_max_eigenvalue(this)
965 type(derivatives_t), intent(in) :: this
966
967 integer :: i
968
969 push_sub(derivatives_lapl_get_max_eigenvalue)
970
971 derivatives_lapl_get_max_eigenvalue = m_zero
972 if ((this%stencil_type == der_star .or. this%stencil_type == der_stargeneral) &
973 .and. .not. this%mesh%use_curvilinear) then
974 ! use Fourier transform of stencil evaluated at the maximum phase
975 do i = 1, this%lapl%stencil%size
976 derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
977 (-1)**maxval(abs(this%lapl%stencil%points(:, i)))*this%lapl%w(i, 1)
978 end do
979 derivatives_lapl_get_max_eigenvalue = abs(derivatives_lapl_get_max_eigenvalue)
980 else
981 ! use upper bound from continuum for other stencils
982 do i = 1, this%dim
983 derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
984 m_pi**2/this%mesh%spacing(i)**2
985 end do
986 end if
988 pop_sub(derivatives_lapl_get_max_eigenvalue)
990
991
992#include "undef.F90"
993#include "real.F90"
994#include "derivatives_inc.F90"
995
996#include "undef.F90"
997#include "complex.F90"
998#include "derivatives_inc.F90"
1000end module derivatives_oct_m
1001
1002!! Local Variables:
1003!! mode: f90
1004!! coding: utf-8
1005!! End:
double sqrt(double __x) __attribute__((__nothrow__
subroutine, public accel_kernel_build(this, file_name, kernel_name, flags)
Definition: accel.F90:2057
pure logical function, public accel_is_enabled()
Definition: accel.F90:403
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:118
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 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, public dderivatives_batch_grad(der, ffb, opffb, ghost_update, set_bc, to_cartesian, factor)
apply the gradient to a batch of mesh functions
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 get_rhs_lapl(coord_system, polynomials, chi, rhs, stencil_primitive_coordinates)
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, public zderivatives_batch_grad(der, ffb, opffb, ghost_update, set_bc, to_cartesian, factor)
apply the gradient to a batch of mesh functions
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 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)