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, chi, rhs, stencil_primitive_coordinates)
641 integer,
intent(in) :: polynomials(:,:)
642 real(real64),
intent(in) :: chi(:)
643 real(real64),
intent(out) :: rhs(:)
644 logical,
intent(in) :: stencil_primitive_coordinates
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)
653 dim =
size(polynomials, dim=1)
655 if (stencil_primitive_coordinates)
then
657 metric_inverse = coord_system%metric_inverse(chi)
658 hessian_trace = coord_system%trace_hessian(chi)
663 metric_inverse(i, i) =
m_one
670 do j = 1,
size(polynomials, dim=2)
674 if (polynomials(i, j) <= 2)
then
675 powers(polynomials(i, j)) = powers(polynomials(i, j)) + 1
681 if (powers(2) == 1 .and. powers(0) == dim - 1)
then
683 if (polynomials(i, j) == 2)
then
684 rhs(j) =
m_two*metric_inverse(i, i)
690 if (powers(1) == 2 .and. powers(0) == dim - 2)
then
692 if (polynomials(i, j) == 1)
then
694 if (polynomials(k, j) == 1)
then
695 rhs(j) = metric_inverse(i, k) + metric_inverse(k, i)
704 if (powers(1) == 1 .and. powers(0) == dim - 1)
then
706 if (polynomials(i, j) == 1)
then
707 rhs(j) = hessian_trace(i)
718 integer,
intent(in) :: polynomials(:,:)
719 integer,
intent(in) :: dir
720 real(real64),
intent(out) :: rhs(:)
727 dim =
size(polynomials, dim=1)
731 do j = 1,
size(polynomials, dim=2)
734 if (k == dir .and. polynomials(k, j) /= 1) this_one = .false.
735 if (k /= dir .and. polynomials(k, j) /= 0) this_one = .false.
737 if (this_one) rhs(j) =
m_one
748 integer,
intent(in) :: nderiv
749 integer,
intent(in) :: ideriv
751 character(len=80),
optional,
intent(in) :: name
752 logical,
optional,
intent(in) :: verbose
753 integer,
optional,
intent(in) :: order
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_
766 safe_allocate(rhs(1:op(1)%stencil%size, nderiv))
767 safe_allocate(pol(1:der%dim, 1:op(1)%stencil%size))
769 if (nderiv == 1)
then
770 if (ideriv <= der%dim)
then
774 call stencil_pol_lapl(der%stencil_type, op(1)%stencil, der%dim, order_, pol)
777 else if (nderiv == der%dim)
then
784 message(1) =
"Error: derivatives_make_discretization can only be called with nderiv = 1 or der%dim."
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))
793 message(1) =
'Info: Generating weights for finite-difference discretization of ' // trim(name_)
798 pow_max = maxval(pol)
799 safe_allocate(powers(1:der%dim, 0:pow_max))
804 if (op(1)%const_w) p_max = 1
809 if (nderiv == 1)
then
810 if (ideriv <= der%dim)
then
814 call get_rhs_lapl(der%mesh%coord_system, pol, der%mesh%chi(:, p), rhs(:,1), der%stencil_primitive_coordinates)
816 else if (nderiv == der%dim)
then
827 do i = 1, op(1)%stencil%size
828 if (ideriv <= der%dim)
then
830 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
833 if (der%stencil_primitive_coordinates)
then
834 x = real(op(1)%stencil%points(:, i), real64)*der%mesh%spacing
836 x = der%mesh%x(p + op(1)%ri(i, op(1)%rimap(p)), :) - der%mesh%x(p, :)
841 x = x*
sqrt(der%masses)
846 powers(:, k) = x*powers(:, k-1)
851 do j = 2, op(1)%stencil%size
852 mat(j, i) = powers(1, pol(1, j))
854 mat(j, i) = mat(j, i)*powers(k, pol(k, j))
866 op(i)%w(:, p) = sol(:, i)
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)
899 overlap = this%comm_method /= blocking
910 class(
space_t),
intent(in) :: space
911 character(len=80),
intent(in) :: name
912 integer,
intent(in) :: order
917 call nl_operator_build(space, this%mesh, op(1), this%mesh%np, const_w = .not. this%mesh%use_curvilinear)
939 logical :: mask(1:this%mesh%np)
940 integer :: ip, is, index
944 do ip = 1, this%mesh%np
946 do is = 1, this%lapl%stencil%size
950 if (index > this%mesh%np + this%mesh%pv%np_ghost)
then
964 real(real64) function derivatives_lapl_get_max_eigenvalue(this)
969 push_sub(derivatives_lapl_get_max_eigenvalue)
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
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)
979 derivatives_lapl_get_max_eigenvalue = abs(derivatives_lapl_get_max_eigenvalue)
983 derivatives_lapl_get_max_eigenvalue = derivatives_lapl_get_max_eigenvalue + &
984 m_pi**2/this%mesh%spacing(i)**2
988 pop_sub(derivatives_lapl_get_max_eigenvalue)
994#include "derivatives_inc.F90"
997#include "complex.F90"
998#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)...
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
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.