42 use,
intrinsic :: iso_fortran_env
117 integer,
parameter :: &
118 DER_BC_ZERO_F = 0, & !< function is zero at the boundaries
122 integer,
parameter,
public :: &
129 integer,
parameter :: &
137 type(boundaries_t) :: boundaries
138 type(mesh_t),
pointer :: mesh => null()
140 integer,
private :: periodic_dim = 0
142 integer :: stencil_type = 0
145 logical,
private :: stencil_primitive_coordinates
147 real(real64),
allocatable :: masses(:)
151 real(real64),
private :: lapl_cutoff =
m_zero
153 type(nl_operator_t),
allocatable,
private :: op(:)
155 type(nl_operator_t),
pointer :: lapl => null()
156 type(nl_operator_t),
pointer :: grad(:) => null()
158 integer,
allocatable :: n_ghost(:)
160 integer,
private :: comm_method = 0
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()
172 type(par_vec_handle_batch_t) :: pv_h
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
183 type(accel_kernel_t) :: kernel_uvw_xyz, kernel_dcurl, kernel_zcurl
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
196 integer :: default_stencil
197 character(len=40) :: flags
198 logical :: stencil_primitive_coordinates
203 assert(.not. coord_system%local_basis .or. coord_system%orthogonal)
207 der%periodic_dim = space%periodic_dim
232 default_stencil = der_star
236 call parse_variable(namespace,
'DerivativesStencil', default_stencil, der%stencil_type)
269 if (
present(order))
then
311 stencil_primitive_coordinates = .not. coord_system%local_basis
312 call parse_variable(namespace,
'StencilPrimitiveCoordinates', stencil_primitive_coordinates, der%stencil_primitive_coordinates)
315 safe_allocate(der%masses(1:space%dim))
319 safe_allocate(der%op(1:der%dim + 1))
321 der%lapl => der%op(der%dim + 1)
327 safe_allocate(der%n_ghost(1:der%dim))
330 der%n_ghost(idir) = maxval(abs(der%lapl%stencil%points(idir, :)))
335 nullify(der%to_coarser)
336 nullify(der%to_finer)
340 select type (coord_system)
344 if (der%dim > 3)
then
345 message(1) =
"Calculation of derivatives on the GPU with dimension > 3 are only implemented for Cartesian coordinates."
348 write(flags,
'(A,I1.1)')
' -DDIMENSION=', der%dim
368 assert(
allocated(der%op))
370 do idim = 1, der%dim+1
374 safe_deallocate_a(der%masses)
376 safe_deallocate_a(der%n_ghost)
378 safe_deallocate_a(der%op)
379 nullify(der%lapl, der%grad)
383 nullify(der%to_coarser)
384 nullify(der%to_finer)
396 class(
space_t),
intent(in) :: space
398 character(len=80),
optional,
intent(in) :: name
399 integer,
optional,
intent(in) :: order
401 character(len=80) :: name_
413 select case (der%stencil_type)
433 real(real64),
intent(out) :: lapl(:)
437 assert(ubound(lapl, dim=1) >= der%mesh%np)
453 character :: dir_label
457 assert(
associated(der%grad))
467 select case (der%stencil_type)
530 subroutine derivatives_build(der, namespace, space, mesh, qvector, regenerate, verbose)
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
542 integer :: np_zero_bc
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))
559 if (mesh%use_curvilinear) const_w_ = .false.
566 safe_deallocate_a(der%op(i)%w)
568 call nl_operator_build(space, mesh, der%op(i), der%mesh%np, const_w = const_w_, &
569 regenerate=regenerate)
573 assert(np_zero_bc > mesh%np .and. np_zero_bc <= mesh%np_part)
575 select case (der%stencil_type)
578 do i = 1, der%dim + 1
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(:, :)
607 select case (stencil_type)
620 integer,
intent(in) :: stencil_type
622 integer,
intent(in) :: dim
623 integer,
intent(in) :: order
624 integer,
intent(inout) :: polynomials(:, :)
626 select case (stencil_type)
639 subroutine get_rhs_lapl(coord_system, polynomials, x, rhs, stencil_primitive_coordinates)
641 integer,
intent(in) :: polynomials(:,:)
642 real(real64),
intent(in) :: x(:)
643 real(real64),
intent(out) :: rhs(:)
644 logical,
intent(in) :: stencil_primitive_coordinates
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)
654 dim =
size(polynomials, dim=1)
656 if (stencil_primitive_coordinates)
then
658 select type (coord => coord_system)
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)))
665 metric(1:dim, 1:dim) = matmul(jacobian(1:dim, 1:dim), transpose(jacobian(1:dim, 1:dim)))
668 message(1) =
"Weight computation not implemented for the coordinate system chosen."
682 do j = 1,
size(polynomials, dim=2)
686 if (polynomials(i, j) <= 2)
then
687 powers(polynomials(i, j)) = powers(polynomials(i, j)) + 1
693 if (powers(2) == 1 .and. powers(0) == dim - 1)
then
695 if (polynomials(i, j) == 2)
then
696 rhs(j) =
m_two*metric(i, i)
702 if (powers(1) == 2 .and. powers(0) == dim - 2)
then
704 if (polynomials(i, j) == 1)
then
706 if (polynomials(k, j) == 1)
then
707 rhs(j) = metric(i, k) + metric(k, i)
716 if (powers(1) == 1 .and. powers(0) == dim - 1)
then
718 if (polynomials(i, j) == 1)
then
719 rhs(j) = hessian_trace(i)
730 integer,
intent(in) :: polynomials(:,:)
731 integer,
intent(in) :: dir
732 real(real64),
intent(out) :: rhs(:)
739 dim =
size(polynomials, dim=1)
743 do j = 1,
size(polynomials, dim=2)
746 if (k == dir .and. polynomials(k, j) /= 1) this_one = .false.
747 if (k /= dir .and. polynomials(k, j) /= 0) this_one = .false.
749 if (this_one) rhs(j) =
m_one
760 integer,
intent(in) :: nderiv
761 integer,
intent(in) :: ideriv
763 character(len=80),
optional,
intent(in) :: name
764 logical,
optional,
intent(in) :: verbose
765 integer,
optional,
intent(in) :: order
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_
778 safe_allocate(rhs(1:op(1)%stencil%size, nderiv))
779 safe_allocate(pol(1:der%dim, 1:op(1)%stencil%size))
781 if (nderiv == 1)
then
782 if (ideriv <= der%dim)
then
786 call stencil_pol_lapl(der%stencil_type, op(1)%stencil, der%dim, order_, pol)
789 else if (nderiv == der%dim)
then
796 message(1) =
"Error: derivatives_make_discretization can only be called with nderiv = 1 or der%dim."
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))
805 message(1) =
'Info: Generating weights for finite-difference discretization of ' // trim(name_)
810 pow_max = maxval(pol)
811 safe_allocate(powers(1:der%dim, 0:pow_max))
816 if (op(1)%const_w) p_max = 1
822 if (nderiv == 1)
then
823 if (ideriv <= der%dim)
then
827 call get_rhs_lapl(der%mesh%coord_system, pol, x, rhs(:,1), der%stencil_primitive_coordinates)
829 else if (nderiv == der%dim)
then
840 do i = 1, op(1)%stencil%size
841 if (ideriv <= der%dim)
then
843 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
846 if (der%stencil_primitive_coordinates)
then
847 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
849 x = der%mesh%x(p + op(1)%ri(i, op(1)%rimap(p)), :) - der%mesh%x(p, :)
854 x = x*
sqrt(der%masses)
859 powers(:, k) = x*powers(:, k-1)
864 do j = 2, op(1)%stencil%size
865 mat(j, i) = powers(1, pol(1, j))
867 mat(j, i) = mat(j, i)*powers(k, pol(k, j))
879 op(i)%w(:, p) = sol(:, i)
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)
912 overlap = this%comm_method /= blocking
923 class(
space_t),
intent(in) :: space
924 character(len=80),
intent(in) :: name
925 integer,
intent(in) :: order
930 call nl_operator_build(space, this%mesh, op(1), this%mesh%np, const_w = .not. this%mesh%use_curvilinear)
952 logical :: mask(1:this%mesh%np)
953 integer :: ip, is, index
957 do ip = 1, this%mesh%np
959 do is = 1, this%lapl%stencil%size
963 if (index > this%mesh%np + this%mesh%pv%np_ghost)
then
977 real(real64) function derivatives_lapl_get_max_eigenvalue(this)
982 push_sub(derivatives_lapl_get_max_eigenvalue)
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
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)
992 derivatives_lapl_get_max_eigenvalue = abs(derivatives_lapl_get_max_eigenvalue)
996 derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
997 m_pi**2/this%mesh%spacing(i)**2
1001 pop_sub(derivatives_lapl_get_max_eigenvalue)
1007#include "derivatives_inc.F90"
1010#include "complex.F90"
1011#include "derivatives_inc.F90"
double sqrt(double __x) __attribute__((__nothrow__
subroutine, public accel_kernel_build(this, file_name, kernel_name, flags)
pure logical function, public accel_is_enabled()
This module implements batches of mesh functions.
This module implements common operations on batches of mesh functions.
Module implementing boundary conditions in Octopus.
subroutine, public boundaries_end(this)
subroutine, public boundaries_init(this, namespace, space, mesh, qvector)
initialize the boundary contitions
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)...
pure subroutine, public curv_gygi_jacobian(this, xx, chi, jac, natoms)
pure subroutine, public curv_gygi_hessian_trace(this, xx, hessian_trace, natoms)
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
real(real64), parameter, public m_zero
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public m_one
This module is intended to contain "only mathematical" functions and procedures.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_obsolete_variable(namespace, name, rep)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
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:
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.
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.
character pure function, public index2axis(idir)
abstract class to describe coordinate systems
handle to transfer data from the start() to finish() calls.
class representing derivatives
Describes mesh distribution to nodes.
data type for non local operators
The class representing the stencil, which is used for non-local mesh operations.