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 21module derivatives_oct_m 22 use accel_oct_m 23 use batch_oct_m 24 use batch_ops_oct_m 25 use boundaries_oct_m 26 use global_oct_m 27 use iso_c_binding 28 use lalg_adv_oct_m 29 use loct_oct_m 30 use math_oct_m 31 use mesh_oct_m 32 use mesh_function_oct_m 33 use messages_oct_m 34 use namespace_oct_m 35 use nl_operator_oct_m 36 use par_vec_oct_m 37 use parser_oct_m 38 use profiling_oct_m 39 use simul_box_oct_m 40 use stencil_cube_oct_m 41 use stencil_star_oct_m 42 use stencil_starplus_oct_m 43 use stencil_stargeneral_oct_m 44 use stencil_variational_oct_m 45 use transfer_table_oct_m 46 use types_oct_m 47 use utils_oct_m 48 use varinfo_oct_m 49 50! debug purposes 51! use io_binary_oct_m 52! use io_function_oct_m 53! use io_oct_m 54! use unit_oct_m 55! use unit_system_oct_m 56 57 implicit none 58 59 private 60 public :: & 61 derivatives_t, & 62 derivatives_nullify, & 63 derivatives_init, & 64 derivatives_end, & 65 derivatives_build, & 66 derivatives_handle_batch_t, & 67 dderivatives_test, & 68 zderivatives_test, & 69 dderivatives_batch_start, & 70 zderivatives_batch_start, & 71 dderivatives_batch_finish, & 72 zderivatives_batch_finish, & 73 dderivatives_batch_perform, & 74 zderivatives_batch_perform, & 75 dderivatives_perform, & 76 zderivatives_perform, & 77 dderivatives_lapl, & 78 zderivatives_lapl, & 79 derivatives_lapl_diag, & 80 dderivatives_grad, & 81 zderivatives_grad, & 82 dderivatives_batch_grad, & 83 zderivatives_batch_grad, & 84 dderivatives_div, & 85 zderivatives_div, & 86 dderivatives_curl, & 87 zderivatives_curl, & 88 dderivatives_batch_curl, & 89 zderivatives_batch_curl, & 90 dderivatives_partial, & 91 zderivatives_partial, & 92 derivatives_get_lapl 93 94 95 integer, parameter :: & 96 DER_BC_ZERO_F = 0, & !< function is zero at the boundaries 97 DER_BC_ZERO_DF = 1, & !< first derivative of the function is zero 98 DER_BC_PERIOD = 2 !< boundary is periodic 99 100 integer, parameter :: & 101 DER_STAR = 1, & 102 DER_VARIATIONAL = 2, & 103 DER_CUBE = 3, & 104 DER_STARPLUS = 4, & 105 DER_STARGENERAL = 5 106 107 integer, parameter :: & 108 BLOCKING = 1, & 109 NON_BLOCKING = 2 110 111 type derivatives_t 112 ! Components are public by default 113 type(boundaries_t) :: boundaries 114 type(mesh_t), pointer :: mesh !< pointer to the underlying mesh 115 integer :: dim !< dimensionality of the space (sb%dim) 116 integer :: order !< order of the discretization (value depends on stencil) 117 integer :: stencil_type !< type of discretization 118 119 FLOAT :: masses(MAX_DIM) !< we can have different weights (masses) per space direction 120 121 !> If the so-called variational discretization is used, this controls a 122 !! possible filter on the Laplacian. 123 FLOAT, private :: lapl_cutoff 124 125 type(nl_operator_t), pointer, private :: op(:) !< op(1:conf%dim) => gradient 126 !! op(conf%dim+1) => Laplacian 127 type(nl_operator_t), pointer :: lapl !< these are just shortcuts for op 128 type(nl_operator_t), pointer :: grad(:) 129 130 integer :: n_ghost(MAX_DIM) !< ghost points to add in each dimension 131#if defined(HAVE_MPI) 132 integer, private :: comm_method 133#endif 134 type(derivatives_t), pointer :: finer 135 type(derivatives_t), pointer :: coarser 136 type(transfer_table_t), pointer :: to_finer 137 type(transfer_table_t), pointer :: to_coarser 138 end type derivatives_t 139 140 type derivatives_handle_batch_t 141 private 142#ifdef HAVE_MPI 143 type(pv_handle_batch_t) :: pv_h 144#endif 145 type(derivatives_t), pointer :: der 146 type(nl_operator_t), pointer :: op 147 type(batch_t), pointer :: ff 148 type(batch_t), pointer :: opff 149 logical :: ghost_update 150 logical :: factor_present 151 FLOAT :: factor 152 end type derivatives_handle_batch_t 153 154 type(accel_kernel_t) :: kernel_uvw_xyz, kernel_dcurl, kernel_zcurl 155 156 type(profile_t), save :: gradient_prof, divergence_prof, curl_prof, batch_gradient_prof, curl_batch_prof 157 158contains 159 160 ! --------------------------------------------------------- 161 elemental subroutine derivatives_nullify(this) 162 type(derivatives_t), intent(out) :: this 163 164 call boundaries_nullify(this%boundaries) 165 nullify(this%mesh) 166 this%dim = 0 167 this%order = 0 168 this%stencil_type = 0 169 this%masses = M_ZERO 170 this%lapl_cutoff = M_ZERO 171 nullify(this%op, this%lapl, this%grad) 172 this%n_ghost = 0 173#if defined(HAVE_MPI) 174 this%comm_method = 0 175#endif 176 nullify(this%finer, this%coarser) 177 nullify(this%to_finer, this%to_coarser) 178 179 end subroutine derivatives_nullify 180 181 ! --------------------------------------------------------- 182 subroutine derivatives_init(der, namespace, sb, use_curvilinear, order) 183 type(derivatives_t), target, intent(out) :: der 184 type(namespace_t), intent(in) :: namespace 185 type(simul_box_t), intent(in) :: sb 186 logical, intent(in) :: use_curvilinear 187 integer, optional, intent(in) :: order 188 189 integer :: idir 190 integer :: default_stencil 191 character(len=40) :: flags 192 193 PUSH_SUB(derivatives_init) 194 195 ! copy this value to my structure 196 der%dim = sb%dim 197 198 !%Variable DerivativesStencil 199 !%Type integer 200 !%Default stencil_star 201 !%Section Mesh::Derivatives 202 !%Description 203 !% Decides what kind of stencil is used, <i>i.e.</i> which points, around 204 !% each point in the mesh, are the neighboring points used in the 205 !% expression of the differential operator. 206 !% 207 !% If curvilinear coordinates are to be used, then only the <tt>stencil_starplus</tt> 208 !% or the <tt>stencil_cube</tt> may be used. We only recommend the <tt>stencil_starplus</tt>, 209 !% since the cube typically needs far too much memory. 210 !%Option stencil_star 1 211 !% A star around each point (<i>i.e.</i>, only points on the axis). 212 !%Option stencil_variational 2 213 !% Same as the star, but with coefficients built in a different way. 214 !%Option stencil_cube 3 215 !% A cube of points around each point. 216 !%Option stencil_starplus 4 217 !% The star, plus a number of off-axis points. 218 !%Option stencil_stargeneral 5 219 !% The general star. Default for non-orthogonal grids. 220 !%End 221 default_stencil = DER_STAR 222 if(use_curvilinear) default_stencil = DER_STARPLUS 223 if(sb%nonorthogonal) default_stencil = DER_STARGENERAL 224 225 call parse_variable(namespace, 'DerivativesStencil', default_stencil, der%stencil_type) 226 227 if(.not.varinfo_valid_option('DerivativesStencil', der%stencil_type)) then 228 call messages_input_error(namespace, 'DerivativesStencil') 229 end if 230 call messages_print_var_option(stdout, "DerivativesStencil", der%stencil_type) 231 232 if(use_curvilinear .and. der%stencil_type < DER_CUBE) call messages_input_error(namespace, 'DerivativesStencil') 233 if(der%stencil_type == DER_VARIATIONAL) then 234 call parse_variable(namespace, 'DerivativesLaplacianFilter', M_ONE, der%lapl_cutoff) 235 end if 236 237 !%Variable DerivativesOrder 238 !%Type integer 239 !%Default 4 240 !%Section Mesh::Derivatives 241 !%Description 242 !% This variable gives the discretization order for the approximation of 243 !% the differential operators. This means, basically, that 244 !% <tt>DerivativesOrder</tt> points are used in each positive/negative 245 !% spatial direction, <i>e.g.</i> <tt>DerivativesOrder = 1</tt> would give 246 !% the well-known three-point formula in 1D. 247 !% The number of points actually used for the Laplacian 248 !% depends on the stencil used. Let <math>O</math> = <tt>DerivativesOrder</tt>, and <math>d</math> = <tt>Dimensions</tt>. 249 !% <ul> 250 !% <li> <tt>stencil_star</tt>: <math>2 O d + 1</math> 251 !% <li> <tt>stencil_cube</tt>: <math>(2 O + 1)^d</math> 252 !% <li> <tt>stencil_starplus</tt>: <math>2 O d + 1 + n</math> with <i>n</i> being 8 253 !% in 2D and 24 in 3D. 254 !% </ul> 255 !%End 256 call parse_variable(namespace, 'DerivativesOrder', 4, der%order) 257 ! overwrite order if given as argument 258 if(present(order)) then 259 der%order = order 260 end if 261 262#ifdef HAVE_MPI 263 !%Variable ParallelizationOfDerivatives 264 !%Type integer 265 !%Default non_blocking 266 !%Section Execution::Parallelization 267 !%Description 268 !% This option selects how the communication of mesh boundaries is performed. 269 !%Option blocking 1 270 !% Blocking communication. 271 !%Option non_blocking 2 272 !% Communication is based on non-blocking point-to-point communication. 273 !%End 274 275 call parse_variable(namespace, 'ParallelizationOfDerivatives', NON_BLOCKING, der%comm_method) 276 277 if(.not. varinfo_valid_option('ParallelizationOfDerivatives', der%comm_method)) then 278 call messages_input_error(namespace, 'ParallelizationOfDerivatives') 279 end if 280 281 call messages_obsolete_variable(namespace, 'OverlapDerivatives', 'ParallelizationOfDerivatives') 282#endif 283 284 ! if needed, der%masses should be initialized in modelmb_particles_init 285 der%masses = M_ONE 286 287 ! construct lapl and grad structures 288 SAFE_ALLOCATE(der%op(1:der%dim + 1)) 289 der%grad => der%op 290 der%lapl => der%op(der%dim + 1) 291 292 call derivatives_get_stencil_lapl(der, sb) 293 call derivatives_get_stencil_grad(der) 294 295 ! find out how many ghost points we need in each dimension 296 der%n_ghost(:) = 0 297 do idir = 1, der%dim 298 der%n_ghost(idir) = maxval(abs(der%lapl%stencil%points(idir, :))) 299 end do 300 301 nullify(der%coarser) 302 nullify(der%finer) 303 nullify(der%to_coarser) 304 nullify(der%to_finer) 305 306 if(accel_is_enabled()) then 307 write(flags, '(A,I1.1)') ' -DDIMENSION=', der%dim 308 call accel_kernel_build(kernel_uvw_xyz, 'uvw_to_xyz.cl', 'uvw_to_xyz', flags) 309 call accel_kernel_build(kernel_dcurl, 'curl.cl', 'dcurl', flags = '-DRTYPE_DOUBLE') 310 call accel_kernel_build(kernel_zcurl, 'curl.cl', 'zcurl', flags = '-DRTYPE_COMPLEX') 311 end if 312 313 POP_SUB(derivatives_init) 314 end subroutine derivatives_init 315 316 317 ! --------------------------------------------------------- 318 subroutine derivatives_end(der) 319 type(derivatives_t), intent(inout) :: der 320 321 integer :: idim 322 323 PUSH_SUB(derivatives_end) 324 325 ASSERT(associated(der%op)) 326 327 do idim = 1, der%dim+1 328 call nl_operator_end(der%op(idim)) 329 end do 330 331 SAFE_DEALLOCATE_P(der%op) 332 nullify(der%lapl, der%grad) 333 334 nullify(der%coarser) 335 nullify(der%finer) 336 nullify(der%to_coarser) 337 nullify(der%to_finer) 338 339 call boundaries_end(der%boundaries) 340 341 POP_SUB(derivatives_end) 342 end subroutine derivatives_end 343 344 345 ! --------------------------------------------------------- 346 subroutine derivatives_get_stencil_lapl(der, sb) 347 type(derivatives_t), intent(inout) :: der 348 type(simul_box_t), optional,intent(in) :: sb 349 350 PUSH_SUB(derivatives_get_stencil_lapl) 351 352 ASSERT(associated(der%lapl)) 353 354 ! initialize nl operator 355 call nl_operator_init(der%lapl, "Laplacian") 356 357 ! create stencil 358 select case(der%stencil_type) 359 case(DER_STAR, DER_VARIATIONAL) 360 call stencil_star_get_lapl(der%lapl%stencil, der%dim, der%order) 361 case(DER_CUBE) 362 call stencil_cube_get_lapl(der%lapl%stencil, der%dim, der%order) 363 case(DER_STARPLUS) 364 call stencil_starplus_get_lapl(der%lapl%stencil, der%dim, der%order) 365 case(DER_STARGENERAL) 366 call stencil_stargeneral_get_arms(der%lapl%stencil, sb) 367 call stencil_stargeneral_get_lapl(der%lapl%stencil, der%dim, der%order) 368 end select 369 370 POP_SUB(derivatives_get_stencil_lapl) 371 372 end subroutine derivatives_get_stencil_lapl 373 374 375 ! --------------------------------------------------------- 376 !> Returns the diagonal elements of the Laplacian, needed for preconditioning 377 subroutine derivatives_lapl_diag(der, lapl) 378 type(derivatives_t), intent(in) :: der 379 FLOAT, intent(out) :: lapl(:) !< lapl(mesh%np) 380 381 PUSH_SUB(derivatives_lapl_diag) 382 383 ASSERT(ubound(lapl, DIM=1) >= der%mesh%np) 384 385 ! the Laplacian is a real operator 386 call dnl_operator_operate_diag(der%lapl, lapl) 387 388 POP_SUB(derivatives_lapl_diag) 389 390 end subroutine derivatives_lapl_diag 391 392 393 ! --------------------------------------------------------- 394 subroutine derivatives_get_stencil_grad(der) 395 type(derivatives_t), intent(inout) :: der 396 397 398 integer :: ii 399 character :: dir_label 400 401 PUSH_SUB(derivatives_get_stencil_grad) 402 403 ASSERT(associated(der%grad)) 404 405 ! initialize nl operator 406 do ii = 1, der%dim 407 dir_label = ' ' 408 if(ii < 5) dir_label = index2axis(ii) 409 410 call nl_operator_init(der%grad(ii), "Gradient "//dir_label) 411 412 ! create stencil 413 select case(der%stencil_type) 414 case(DER_STAR, DER_VARIATIONAL) 415 call stencil_star_get_grad(der%grad(ii)%stencil, ii, der%order) 416 case(DER_CUBE) 417 call stencil_cube_get_grad(der%grad(ii)%stencil, der%dim, der%order) 418 case(DER_STARPLUS) 419 call stencil_starplus_get_grad(der%grad(ii)%stencil, der%dim, ii, der%order) 420 case(DER_STARGENERAL) 421 ! use the simple star stencil 422 call stencil_star_get_grad(der%grad(ii)%stencil, ii, der%order) 423 end select 424 end do 425 426 POP_SUB(derivatives_get_stencil_grad) 427 428 end subroutine derivatives_get_stencil_grad 429 430 ! --------------------------------------------------------- 431 subroutine derivatives_build(der, namespace, mesh) 432 type(derivatives_t), intent(inout) :: der 433 type(namespace_t), intent(in) :: namespace 434 type(mesh_t), target, intent(in) :: mesh 435 436 integer, allocatable :: polynomials(:,:) 437 FLOAT, allocatable :: rhs(:,:) 438 integer :: i 439 logical :: const_w_ 440 character(len=32) :: name 441 type(nl_operator_t) :: auxop 442 integer :: np_zero_bc 443 444 PUSH_SUB(derivatives_build) 445 446 call boundaries_init(der%boundaries, namespace, mesh) 447 448 ASSERT(associated(der%op)) 449 ASSERT(der%stencil_type>=DER_STAR .and. der%stencil_type<=DER_STARGENERAL) 450 ASSERT(.not.(der%stencil_type==DER_VARIATIONAL .and. mesh%use_curvilinear)) 451 452 der%mesh => mesh ! make a pointer to the underlying mesh 453 454 const_w_ = .true. 455 456 ! need non-constant weights for curvilinear and scattering meshes 457 if(mesh%use_curvilinear) const_w_ = .false. 458 459 np_zero_bc = 0 460 461 ! build operators 462 do i = 1, der%dim+1 463 call nl_operator_build(mesh, der%op(i), der%mesh%np, const_w = const_w_) 464 np_zero_bc = max(np_zero_bc, nl_operator_np_zero_bc(der%op(i))) 465 end do 466 467 ASSERT(np_zero_bc > mesh%np .and. np_zero_bc <= mesh%np_part) 468 469 select case(der%stencil_type) 470 471 case(DER_STAR) ! Laplacian and gradient have different stencils 472 do i = 1, der%dim + 1 473 SAFE_ALLOCATE(polynomials(1:der%dim, 1:der%op(i)%stencil%npoly)) 474 SAFE_ALLOCATE(rhs(1:der%op(i)%stencil%size, 1:1)) 475 476 if(i <= der%dim) then ! gradient 477 call stencil_star_polynomials_grad(i, der%order, polynomials) 478 call get_rhs_grad(i, rhs(:,1)) 479 name = index2axis(i) // "-gradient" 480 else ! Laplacian 481 call stencil_star_polynomials_lapl(der%dim, der%order, polynomials) 482 call get_rhs_lapl(rhs(:,1)) 483 name = "Laplacian" 484 end if 485 486 call derivatives_make_discretization(der%dim, der%mesh, der%masses, polynomials, rhs, 1, der%op(i:i), name) 487 SAFE_DEALLOCATE_A(polynomials) 488 SAFE_DEALLOCATE_A(rhs) 489 end do 490 491 case(DER_CUBE) ! Laplacian and gradient have similar stencils 492 493 SAFE_ALLOCATE(polynomials(1:der%dim, 1:der%op(1)%stencil%npoly)) 494 SAFE_ALLOCATE(rhs(1:der%op(1)%stencil%size, 1:der%dim + 1)) 495 call stencil_cube_polynomials_lapl(der%dim, der%order, polynomials) 496 497 do i = 1, der%dim 498 call get_rhs_grad(i, rhs(:,i)) 499 end do 500 call get_rhs_lapl(rhs(:, der%dim+1)) 501 502 name = "derivatives" 503 call derivatives_make_discretization(der%dim, der%mesh, der%masses, polynomials, rhs, der%dim+1, der%op(:), name) 504 505 SAFE_DEALLOCATE_A(polynomials) 506 SAFE_DEALLOCATE_A(rhs) 507 508 case(DER_STARPLUS) 509 do i = 1, der%dim 510 SAFE_ALLOCATE(polynomials(1:der%dim, 1:der%op(i)%stencil%npoly)) 511 SAFE_ALLOCATE(rhs(1:der%op(i)%stencil%size, 1:1)) 512 call stencil_starplus_pol_grad(der%dim, i, der%order, polynomials) 513 call get_rhs_grad(i, rhs(:, 1)) 514 name = index2axis(i) // "-gradient" 515 call derivatives_make_discretization(der%dim, der%mesh, der%masses, polynomials, rhs, 1, der%op(i:i), name) 516 SAFE_DEALLOCATE_A(polynomials) 517 SAFE_DEALLOCATE_A(rhs) 518 end do 519 SAFE_ALLOCATE(polynomials(1:der%dim, 1:der%op(der%dim+1)%stencil%npoly)) 520 SAFE_ALLOCATE(rhs(1:der%op(i)%stencil%size, 1:1)) 521 call stencil_starplus_pol_lapl(der%dim, der%order, polynomials) 522 call get_rhs_lapl(rhs(:, 1)) 523 name = "Laplacian" 524 call derivatives_make_discretization(der%dim, der%mesh, der%masses, polynomials, rhs, 1, der%op(der%dim+1:der%dim+1), name) 525 SAFE_DEALLOCATE_A(polynomials) 526 SAFE_DEALLOCATE_A(rhs) 527 528 case(DER_STARGENERAL) 529 530 do i = 1, der%dim 531 der%op(i)%stencil%npoly = der%op(i)%stencil%size ! for gradients we are fine 532 533 SAFE_ALLOCATE(polynomials(1:der%dim, 1:der%op(i)%stencil%npoly)) 534 SAFE_ALLOCATE(rhs(1:der%op(i)%stencil%size, 1:1)) 535 ! use simple star stencil polynomials 536 call stencil_star_polynomials_grad(i, der%order, polynomials) 537 call get_rhs_grad(i, rhs(:, 1)) 538 name = index2axis(i) // "-gradient" 539 ! For directional derivatives the weights are the same as in the orthogonal case. 540 ! Forcing the orthogonal case avoid incurring in ill-defined cases. 541 ! NOTE: this is not so clean and also am not 100% sure is correct. It has to be tested. UDG 542 call derivatives_make_discretization(der%dim, der%mesh, der%masses, polynomials, rhs, 1,& 543 der%op(i:i), name, force_orthogonal = .true.) 544 SAFE_DEALLOCATE_A(polynomials) 545 SAFE_DEALLOCATE_A(rhs) 546 end do 547 548 der%op(der%dim+1)%stencil%npoly = der%op(der%dim+1)%stencil%size & 549 + der%order*(2*der%order-1)*der%op(der%dim+1)%stencil%stargeneral%narms 550 551 SAFE_ALLOCATE(polynomials(1:der%dim, 1:der%op(der%dim+1)%stencil%npoly)) 552 SAFE_ALLOCATE(rhs(1:der%op(der%dim+1)%stencil%size, 1:1)) 553 call stencil_stargeneral_pol_lapl(der%op(der%dim+1)%stencil, der%dim, der%order, polynomials) 554 call get_rhs_lapl(rhs(:, 1)) 555 name = "Laplacian" 556 call derivatives_make_discretization(der%dim, der%mesh, der%masses, polynomials, rhs, 1, der%op(der%dim+1:der%dim+1), name) 557 SAFE_DEALLOCATE_A(polynomials) 558 SAFE_DEALLOCATE_A(rhs) 559 560 case(DER_VARIATIONAL) 561 ! we have the explicit coefficients 562 call stencil_variational_coeff_lapl(der%dim, der%order, mesh%spacing, der%lapl, alpha = der%lapl_cutoff) 563 564 end select 565 566 567 568 ! Here the Laplacian is forced to be self-adjoint, and the gradient to be skew-self-adjoint 569 if(mesh%use_curvilinear .and. (.not. der%mesh%sb%mr_flag)) then 570 do i = 1, der%dim 571 call nl_operator_init(auxop, "auxop") 572 call nl_operator_skewadjoint(der%grad(i), auxop, der%mesh) 573 574 call nl_operator_end(der%grad(i)) 575 call nl_operator_copy(der%grad(i), auxop) 576 call nl_operator_end(auxop) 577 end do 578 call nl_operator_init(auxop, "auxop") 579 call nl_operator_selfadjoint(der%lapl, auxop, der%mesh) 580 581 call nl_operator_end(der%lapl) 582 call nl_operator_copy(der%lapl, auxop) 583 call nl_operator_end(auxop) 584 end if 585 586 ! useful for debug purposes 587! call dderivatives_test(der, 1, 1, 1) 588! call exit(1) 589 590 591 POP_SUB(derivatives_build) 592 593 contains 594 595 ! --------------------------------------------------------- 596 subroutine get_rhs_lapl(rhs) 597 FLOAT, intent(out) :: rhs(:) 598 599 integer :: i, j, k 600 logical :: this_one 601 602 PUSH_SUB(derivatives_build.get_rhs_lapl) 603 604 ! find right-hand side for operator 605 rhs(:) = M_ZERO 606 do i = 1, der%dim 607 do j = 1, der%lapl%stencil%npoly 608 this_one = .true. 609 do k = 1, der%dim 610 if(k == i .and. polynomials(k, j) /= 2) this_one = .false. 611 if(k /= i .and. polynomials(k, j) /= 0) this_one = .false. 612 end do 613 if(this_one) rhs(j) = M_TWO 614 end do 615 end do 616 617 618 POP_SUB(derivatives_build.get_rhs_lapl) 619 end subroutine get_rhs_lapl 620 621 ! --------------------------------------------------------- 622 subroutine get_rhs_grad(dir, rhs) 623 integer, intent(in) :: dir 624 FLOAT, intent(out) :: rhs(:) 625 626 integer :: j, k 627 logical :: this_one 628 629 PUSH_SUB(derivatives_build.get_rhs_grad) 630 631 ! find right-hand side for operator 632 rhs(:) = M_ZERO 633 do j = 1, der%grad(dir)%stencil%npoly 634 this_one = .true. 635 do k = 1, der%dim 636 if(k == dir .and. polynomials(k, j) /= 1) this_one = .false. 637 if(k /= dir .and. polynomials(k, j) /= 0) this_one = .false. 638 end do 639 if(this_one) rhs(j) = M_ONE 640 end do 641 642 POP_SUB(derivatives_build.get_rhs_grad) 643 end subroutine get_rhs_grad 644 645 end subroutine derivatives_build 646 647 648 ! --------------------------------------------------------- 649 subroutine derivatives_make_discretization(dim, mesh, masses, pol, rhs, nderiv, op, name, force_orthogonal) 650 integer, intent(in) :: dim 651 type(mesh_t), intent(in) :: mesh 652 FLOAT, intent(in) :: masses(:) 653 integer, intent(in) :: pol(:,:) 654 integer, intent(in) :: nderiv 655 FLOAT, intent(inout) :: rhs(:,:) 656 type(nl_operator_t), intent(inout) :: op(:) 657 character(len=32), intent(in) :: name 658 logical, optional, intent(in) :: force_orthogonal 659 660 integer :: p, p_max, i, j, k, pow_max 661 FLOAT :: x(MAX_DIM) 662 FLOAT, allocatable :: mat(:,:), sol(:,:), powers(:,:) 663 664 PUSH_SUB(derivatives_make_discretization) 665 666 SAFE_ALLOCATE(mat(1:op(1)%stencil%npoly, 1:op(1)%stencil%size)) 667 SAFE_ALLOCATE(sol(1:op(1)%stencil%size, 1:nderiv)) 668 669 message(1) = 'Info: Generating weights for finite-difference discretization of ' // trim(name) 670 call messages_info(1) 671 672 ! use to generate power lookup table 673 pow_max = maxval(pol) 674 SAFE_ALLOCATE(powers(1:dim, 0:pow_max)) 675 powers(:,:) = M_ZERO 676 powers(:,0) = M_ONE 677 678 p_max = op(1)%np 679 if(op(1)%const_w) p_max = 1 680 681 do p = 1, p_max 682 ! first polynomial is just a constant 683 mat(1,:) = M_ONE 684 ! i indexes the point in the stencil 685 do i = 1, op(1)%stencil%size 686 if(mesh%use_curvilinear) then 687 do j = 1, dim 688 x(j) = mesh%x(p + op(1)%ri(i, op(1)%rimap(p)), j) - mesh%x(p, j) 689 end do 690 else 691 do j = 1, dim 692 x(j) = TOFLOAT(op(1)%stencil%points(j, i))*mesh%spacing(j) 693 end do 694 ! TODO : this internal if clause is inefficient - the condition is determined globally 695 if (mesh%sb%nonorthogonal .and. .not. optional_default(force_orthogonal, .false.)) & 696 x(1:dim) = matmul(mesh%sb%rlattice_primitive(1:dim,1:dim), x(1:dim)) 697 end if 698 699! NB: these masses are applied on the cartesian directions. Should add a check for non-orthogonal axes 700 do j = 1, dim 701 x(j) = x(j)*sqrt(masses(j)) 702 end do 703 704 ! calculate powers 705 do j = 1, dim 706 powers(j, 1) = x(j) 707 do k = 2, pow_max 708 powers(j, k) = x(j)*powers(j, k-1) 709 end do 710 end do 711 712 ! generate the matrix 713 ! j indexes the polynomial being used 714 do j = 2, op(1)%stencil%npoly 715 mat(j, i) = powers(1, pol(1, j)) 716 do k = 2, dim 717 mat(j, i) = mat(j, i)*powers(k, pol(k, j)) 718 end do 719 end do 720 end do ! loop over i = point in stencil 721 722 ! linear problem to solve for derivative weights: 723 ! mat * sol = rhs 724 725 if (op(1)%stencil%npoly==op(1)%stencil%size) then 726 call lalg_linsyssolve(op(1)%stencil%size, nderiv, mat, rhs, sol) 727 else 728 call lalg_svd_inverse(op(1)%stencil%npoly, op(1)%stencil%size, mat, CNST(1.e-10)) 729 sol = matmul(transpose(mat(1:op(1)%stencil%size, 1:op(1)%stencil%size)), rhs) 730 end if 731 732 do i = 1, nderiv 733!MJV 10 nov 2016: changed n to i below in sol(:,n): was only erroneous in cube case when several 734!derivatives are calculated in 1 batch by this subroutine. All weights contained 735!the last derivative`s weights (gradient z) 736!op(1) is used systematically above to get dimensions, but here we have to save 737!all operator stencil weights 738! once we are happy and convinced, remove this comment 739 op(i)%w(:, p) = sol(:, i) 740 end do 741 742 end do ! loop over points p 743 744 do i = 1, nderiv 745 call nl_operator_update_weights(op(i)) 746 end do 747 748 SAFE_DEALLOCATE_A(mat) 749 SAFE_DEALLOCATE_A(sol) 750 SAFE_DEALLOCATE_A(powers) 751 752 POP_SUB(derivatives_make_discretization) 753 end subroutine derivatives_make_discretization 754 755#ifdef HAVE_MPI 756 ! --------------------------------------------------------- 757 logical function derivatives_overlap(this) result(overlap) 758 type(derivatives_t), intent(in) :: this 759 760 PUSH_SUB(derivatives_overlap) 761 762 overlap = this%comm_method /= BLOCKING 763 764 POP_SUB(derivatives_overlap) 765 end function derivatives_overlap 766#endif 767 768 ! --------------------------------------------------------- 769 subroutine derivatives_get_lapl(this, op, name, order) 770 type(derivatives_t), intent(in) :: this 771 type(nl_operator_t), intent(inout) :: op(:) 772 character(len=32), intent(in) :: name 773 integer, intent(in) :: order 774 775 integer, allocatable :: polynomials(:,:) 776 FLOAT, allocatable :: rhs(:,:) 777 integer :: i, j, k 778 logical :: this_one 779 780 PUSH_SUB(derivatives_get_lapl) 781 782 call nl_operator_init(op(1), name) 783 if(this%mesh%sb%nonorthogonal) then 784 call stencil_stargeneral_get_arms(op(1)%stencil, this%mesh%sb) 785 call stencil_stargeneral_get_lapl(op(1)%stencil, this%dim, order) 786 else 787 call stencil_star_get_lapl(op(1)%stencil, this%dim, order) 788 end if 789 call nl_operator_build(this%mesh, op(1), this%mesh%np, const_w = .not. this%mesh%use_curvilinear) 790 791 !At the moment this code is almost copy-pasted from derivatives_build. 792 if(this%mesh%sb%nonorthogonal) then 793 op(1)%stencil%npoly = op(1)%stencil%size & 794 + order*(2*order-1)*op(1)%stencil%stargeneral%narms 795 end if 796 SAFE_ALLOCATE(polynomials(1:this%dim, 1:op(1)%stencil%npoly)) 797 SAFE_ALLOCATE(rhs(1:op(1)%stencil%size, 1:1)) 798 if(this%mesh%sb%nonorthogonal) then 799 call stencil_stargeneral_pol_lapl(op(1)%stencil, this%dim, order, polynomials) 800 else 801 call stencil_star_polynomials_lapl(this%dim, order, polynomials) 802 end if 803 ! find right-hand side for operator 804 rhs(:,1) = M_ZERO 805 do i = 1, this%dim 806 do j = 1, op(1)%stencil%npoly 807 this_one = .true. 808 do k = 1, this%dim 809 if(k == i .and. polynomials(k, j) /= 2) this_one = .false. 810 if(k /= i .and. polynomials(k, j) /= 0) this_one = .false. 811 end do 812 if(this_one) rhs(j,1) = M_TWO 813 end do 814 end do 815 call derivatives_make_discretization(this%dim, this%mesh, this%masses, & 816 polynomials, rhs, 1, op(1:1), name, force_orthogonal = .true.) 817 SAFE_DEALLOCATE_A(polynomials) 818 SAFE_DEALLOCATE_A(rhs) 819 820 POP_SUB(derivatives_get_lapl) 821 end subroutine derivatives_get_lapl 822 823 824#include "undef.F90" 825#include "real.F90" 826#include "derivatives_inc.F90" 827 828#include "undef.F90" 829#include "complex.F90" 830#include "derivatives_inc.F90" 831 832end module derivatives_oct_m 833 834!! Local Variables: 835!! mode: f90 836!! coding: utf-8 837!! End: 838