Builds a Coriolis block which can be used to solve the inverse problem for the linear components of the plume momentum equations. The result can only be used with fields having the same grid as the template.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=r8), | intent(in) | :: | phi | The dimensionless coriolis parameter |
||
| real(kind=r8), | intent(in) | :: | nu | The dimensionless eddy diffusivity |
||
| integer, | intent(in) | :: | velbound | Location code for the velocity's boundary condition. 1 indicates upper boundary, -1 indicates lower boundary. |
||
| integer, | intent(in) | :: | dvelbound | Location code for the velocity's boundary condition. 1 indicates upper boundary, -1 indicates lower boundary. |
||
| integer, | intent(in) | :: | integrate_bound | Location code for the boundary to perform integrations from. This should be the opposite boundary from where boundary data is stored. |
||
| class(abstract_field), | intent(in) | :: | template | A scalar field with the same grid as any fields passed as arguments to the solve_for method. |
function constructor(phi, nu, velbound, dvelbound, integrate_bound, template) &
result(this)
!* Author: Chris MacMackin
! Date: January 2018
!
! Builds a Coriolis block which can be used to solve the inverse
! problem for the linear components of the plume momentum
! equations. The result can only be used with fields having the
! same grid as the template.
!
real(r8), intent(in) :: phi
!! The dimensionless coriolis parameter
real(r8), intent(in) :: nu
!! The dimensionless eddy diffusivity
integer, intent(in) :: velbound
!! Location code for the velocity's boundary condition. 1
!! indicates upper boundary, -1 indicates lower boundary.
integer, intent(in) :: dvelbound
!! Location code for the velocity's boundary condition. 1
!! indicates upper boundary, -1 indicates lower boundary.
integer, intent(in) :: integrate_bound
!! Location code for the boundary to perform integrations
!! from. This should be the opposite boundary from where
!! boundary data is stored.
class(abstract_field), intent(in) :: template
!! A scalar field with the same grid as any fields passed as
!! arguments to the [[pseudospec_block(type):solve_for]] method.
type(coriolis_block) :: this
type(cheb1d_scalar_field) :: xvals
integer :: i, j, info
real(r8) :: alpha, beta, rcond, s
real(r8), dimension(:,:), allocatable :: domain
complex(r8), dimension(4) :: dummy_in
complex(r8), dimension(4) :: dummy_out
character(len=:), allocatable :: msg
call template%guard_temp()
zero = uniform_scalar_field(0._r8)
this%integrator = pseudospec_block(template)
domain = template%domain()
this%vel_bound_loc = velbound
this%dvel_bound_loc = dvelbound
if (velbound == 1) then
this%xbounds(1:2) = domain(1,2)
else if (velbound == -1) then
this%xbounds(1:2) = domain(1,1)
else
call logger%fatal('coriolis_block', &
'Only boundary location codes 1 or -1 are accepted.')
error stop
end if
if (dvelbound == 1) then
this%xbounds(3:4) = domain(1,2)
else if (dvelbound == -1) then
this%xbounds(3:4) = domain(1,1)
else
call logger%fatal('coriolis_block', &
'Only boundary location codes 1 or -1 are accepted.')
error stop
end if
this%integrate_bound = integrate_bound
xvals = cheb1d_scalar_field(template%elements(), linear, domain(1,1), &
domain(1,2))
! Compute the diagonal elements of \(\bm{B}\)
alpha = sqrt(abs(2*phi/nu))
beta = sqrt(abs(0.5_r8*phi/nu))
s = sign(1._r8,phi)
this%D_r = [-beta, -beta, beta, beta]
this%D_i = [-beta, beta, -beta, beta]
! Use the associations as temporary work arrays, prior to
! computing their actual values
associate (emDx_r => this%eDx_r, emDx_i => this%eDx_i, &
Vinv_r => this%V_r, Vinv_i => this%V_i)
! Compute \(e^{-\bm{D}x}\). Need to do this in two steps to avoid a compiler bug
emDx_i = [(exp(-this%D_r(i)*xvals), i=1,4)]
emDx_r = [(emDx_i(i)*cos(-this%D_i(i)*xvals), i=1,4)]
emDx_i = [(emDx_i(i)*sin(-this%D_i(i)*xvals), i=1,4)]
! Compute \(\bm{V}^{-1}\)
Vinv_r = 0.25_r8*reshape([-s*beta, -s*beta, s*beta, s*beta, &
-beta, -beta, beta, beta, &
0._r8, 0._r8, 0._r8, 0._r8, &
1._r8, 1._r8, 1._r8, 1._r8], &
[4,4])
Vinv_i = 0.25_r8*reshape([s*beta, -s*beta, s*beta, -s*beta, &
-beta, beta, -beta, beta, &
-s, s, s, -s, &
0._r8, 0._r8, 0._r8, 0._r8], &
[4,4])
! Compute \(e^{-\bm{D}x}\bm{V}^{-1}\)
this%emDxVinv_r = reshape([((emDx_r(i)*Vinv_r(i,j) - emDx_i(i)*Vinv_i(i,j), &
i=1,4), j=1,4)], [4, 4])
this%emDxVinv_i = reshape([((emDx_r(i)*Vinv_i(i,j) + emDx_i(i)*Vinv_r(i,j), &
i=1,4), j=1,4)], [4, 4])
end associate
! Compute \(e^{\bm{D}x}\). Need to do this in two steps to avoid a compiler bug
this%eDx_i = [(exp(this%D_r(i)*xvals), i=1,4)]
this%eDx_r = [(this%eDx_i(i)*cos(this%D_i(i)*xvals), i=1,4)]
this%eDx_i = [(this%eDx_i(i)*sin(this%D_i(i)*xvals), i=1,4)]
! Compute \(\bm{V}\)
this%V_r = reshape([-s/alpha, -1._r8/alpha, 0._r8, 1._r8, &
-s/alpha, -1._r8/alpha, 0._r8, 1._r8, &
s/alpha, 1._r8/alpha, 0._r8, 1._r8, &
s/alpha, 1._r8/alpha, 0._r8, 1._r8], &
[4,4])
this%V_i = reshape([-s/alpha, 1._r8/alpha, s, 0._r8, &
s/alpha, -1._r8/alpha, -s, 0._r8, &
-s/alpha, 1._r8/alpha, -s, 0._r8, &
s/alpha, -1._r8/alpha, s, 0._r8], &
[4,4])
! Construct and factor matrix used for satisfying boundary conditions
dummy_in = [(1,0), (0,1), (0.5, 0.5), (-0.5, 0.5)]
this%bound_matrix = reshape([((cmplx(this%V_r(i,j), this%V_i(i,j), r8)* &
exp(cmplx(this%D_r(j), this%D_i(j), r8)* &
this%xbounds(i)), i=1,4), j=1,4)], [4,4])
this%bound_matrix_scaled = this%bound_matrix
call la_gesvx(this%bound_matrix_scaled, dummy_in, dummy_out, &
this%factored_matrix, this%pivots, 'E', trans, &
this%equed, this%r_scales, this%c_scales, &
rcond=rcond, info=info)
if (info /= 0) then
msg = 'Tridiagonal matrix solver returned with flag '//str(info)
call logger%error('coriolis_block',msg)
end if
msg = 'Boundary matrix factored with estimated condition '// &
'number '//str(1._r8/rcond)
call logger%trivia('coriolis_block',msg)
#ifdef DEBUG
call logger%debug('coriolis_block', &
'Successfully constructed Coriolis preconditioner block.')
#endif
call template%clean_temp()
contains
pure function linear(x) result(scalar)
real(r8), dimension(:), intent(in) :: x
!! The position at which this function is evaluated
real(r8) :: scalar
scalar = x(1)
end function linear
end function constructor