Computes the ice shelf velocity at the current time with the current ice thickness.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(ice_shelf), | intent(inout) | :: | this | |||
| class(scalar_field), | intent(in) | :: | basal_drag | A paramter, e.g. coefficient of friction, needed to calculate the drag on basal surface of the glacier. |
||
| logical, | intent(out) | :: | success | True if the integration is successful, false otherwise |
subroutine shelf_solve_velocity(this, basal_drag, success)
!* Author: Chris MacMackin
! Date: May 2017
!
! Computes the ice shelf velocity at the current time with the
! current ice thickness.
!
class(ice_shelf), intent(inout) :: this
class(scalar_field), intent(in) :: basal_drag
!! A paramter, e.g. coefficient of friction, needed to calculate
!! the drag on basal surface of the glacier.
logical, intent(out) :: success
!! True if the integration is successful, false otherwise
integer, save :: nval, kdmax = 20
real(r8), dimension(:), allocatable :: state
integer, dimension(10) :: input
integer, dimension(6) :: info
real(r8), dimension(:), allocatable, save :: work
real(r8), dimension(1) :: real_param
integer, dimension(1) :: int_param
integer :: flag
call basal_drag%guard_temp()
nval = this%velocity%raw_size()
if (allocated(work)) then
if (size(work) < nval*(kdmax+5) + kdmax*(kdmax+3)) then
deallocate(work)
allocate(work(nval*(kdmax+5) + kdmax*(kdmax+3)))
end if
else
allocate(work(nval*(kdmax+5) + kdmax*(kdmax+3)))
end if
input = 0
input(4) = kdmax
input(5) = 1
input(9) = -1
input(10) = 3
state = this%velocity%raw()
call nitsol(nval, state, nitsol_residual, nitsol_precondition, &
1.e-10_r8*nval, 1.e-10_r8*nval, input, info, work, &
real_param, int_param, flag, ddot, dnrm2)
call this%velocity%set_from_raw(state)
this%eta = this%viscosity_law%ice_viscosity(this%velocity, this%ice_temperature(), &
this%time)
this%stale_jacobian = .true.
select case(flag)
case(0)
call logger%trivia('ice_shelf%solve_velocity','Found ice velocity at time '// &
trim(str(this%time)))
success = .true.
case(1)
call logger%error('ice_shelf%solve_velocity','Reached maximum number of'// &
' iterations finding ice velocity')
success = .false.
case default
call logger%error('ice_shelf%solve_velocity','NITSOL failed when finding'// &
' ice velocity with error code '//trim(str(flag)))
success = .false.
end select
call basal_drag%clean_temp()
contains
subroutine nitsol_residual(n, xcur, fcur, rpar, ipar, itrmf)
!! A routine matching the interface expected by NITSOL which
!! returns the residual for the glacier.
integer, intent(in) :: n
!! Dimension of the problem
real(r8), dimension(n), intent(in) :: xcur
!! Array of length `n` containing the current \(x\) value
real(r8), dimension(n), intent(out) :: fcur
!! Array of length `n` containing f(xcur) on output
real(r8), dimension(*), intent(inout) :: rpar
!! Parameter/work array
integer, dimension(*), intent(inout) :: ipar
!! Parameter/work array
integer, intent(out) :: itrmf
!! Termination flag. 0 means normal termination, 1 means
!! failure to produce f(xcur)
type(cheb1d_scalar_field) :: scalar_tmp
integer :: start, finish, bounds_start, bounds_finish
integer, dimension(:), allocatable :: lower, upper
real(r8), dimension(:), allocatable :: bounds
call this%velocity%set_from_raw(xcur)
this%stale_jacobian = .true.
associate(h => this%thickness, uvec => this%velocity, chi => this%chi, &
zeta => this%zeta, eta => this%eta)
eta = this%viscosity_law%ice_viscosity(uvec, this%ice_temperature(), &
this%time)
bounds = this%boundaries%boundary_residuals(h, uvec, eta, this%time)
scalar_tmp = uvec%component(1)
scalar_tmp = 4.0_r8*eta*h*scalar_tmp%d_dx(1)
scalar_tmp = -2.0_r8*chi*h*h%d_dx(1) + scalar_tmp%d_dx(1)
lower = this%boundaries%velocity_lower_bound()
upper = this%boundaries%velocity_upper_bound()
! TODO: Figure out how to make this independent of order which
! values are stored in the field
start = 1
finish = start + this%velocity_upper_bound_size - 1
bounds_start = this%thickness_lower_bound_size &
+ this%thickness_upper_bound_size &
+ this%velocity_lower_bound_size + 1
bounds_finish = bounds_start + this%velocity_upper_bound_size - 1
fcur(start:finish) = bounds(bounds_start:bounds_finish)
start = finish + 1
finish = start + scalar_tmp%raw_size(lower,upper) - 1
fcur(start:finish) = scalar_tmp%raw(lower,upper)
start = finish + 1
finish = start + this%velocity_lower_bound_size - 1
bounds_start = this%thickness_lower_bound_size &
+ this%thickness_upper_bound_size + 1
bounds_finish = bounds_start + this%velocity_lower_bound_size - 1
fcur(start:finish) = bounds(bounds_start:bounds_finish)
end associate
itrmf = 0
end subroutine nitsol_residual
subroutine nitsol_precondition(n, xcur, fcur, ijob, v, z, rpar, ipar, itrmjv)
!! A subroutine matching the interface expected by NITSOL, which
!! acts as a preconditioner.
integer, intent(in) :: n
! Dimension of the problem
real(r8), dimension(n), intent(in) :: xcur
! Array of lenght `n` containing the current $x$ value
real(r8), dimension(n), intent(in) :: fcur
! Array of lenght `n` containing the current \(f(x)\) value
integer, intent(in) :: ijob
! Integer flat indicating which product is desired. 0
! indicates \(z = J\vec{v}\). 1 indicates \(z = P^{-1}\vec{v}\).
real(r8), dimension(n), intent(in) :: v
! An array of length `n` to be multiplied
real(r8), dimension(n), intent(out) :: z
! An array of length n containing the desired product on
! output.
real(r8), dimension(*), intent(inout) :: rpar
! Parameter/work array
integer, dimension(*), intent(inout) :: ipar
! Parameter/work array
integer, intent(out) :: itrmjv
! Termination flag. 0 indcates normal termination, 1
! indicatesfailure to prodce $J\vec{v}$, and 2 indicates
! failure to produce \(P^{-1}\vec{v}\)
type(cheb1d_scalar_field) :: delta_u, tmp
integer :: i, sl, el, su, eu
integer, dimension(2) :: upper_type, lower_type
integer, dimension(:), allocatable :: boundary_types, boundary_locations
real(r8) :: eta_val
if (ijob /= 1) then
itrmjv = 0
return
end if
call delta_u%assign_meta_data(this%velocity)
call delta_u%set_from_raw(v)
sl = this%velocity_size - this%velocity_lower_bound_size + 1
el = this%velocity_size
su = 1
eu = this%velocity_upper_bound_size
boundary_locations = [(i, i=sl,el), (i, i=su,eu)]
upper_type = this%boundaries%velocity_upper_type()
lower_type = this%boundaries%velocity_lower_type()
boundary_types = [(lower_type(1), i=sl,el), (upper_type(1), i=su,eu)]
associate(h => this%thickness, chi => this%chi, zeta => this%zeta, &
jac => this%velocity_jacobian, eta => this%eta)
if (this%stale_jacobian) then
! Mathematically, it should be 4._r8*eta*h which I pass, but
! this sometimes results in an ill-conditioned Jacobian (for
! reasons I'm not clear on). If eta ~ 1 then turns out I get
! good results just ignoring it.
select type(visc => this%viscosity_law)
class is(newtonian_viscosity)
eta_val = eta%get_element(1)
if (eta_val <= 0._r8) eta_val = 1._r8
jac = jacobian_block(4._r8*eta*h, 1, 1, boundary_locs=boundary_locations, &
boundary_types=boundary_types)
class default
jac = jacobian_block(4._r8*h, 1, 1, boundary_locs=boundary_locations, &
boundary_types=boundary_types)
end select
this%stale_jacobian = .false.
end if
delta_u = jac%solve_for(delta_u)
z(1:n) = delta_u%raw()
end associate
itrmjv = 0
end subroutine nitsol_precondition
end subroutine shelf_solve_velocity