Solves the state of the plume for the specified ice properties, at the specified time. This is done using the a quasilinearisation method and a GMRES iterative linear solver.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(plume), | intent(inout) | :: | this | |||
| class(scalar_field), | intent(in) | :: | ice_thickness | Thickness of the ice above the basal surface |
||
| real(kind=r8), | intent(in) | :: | ice_density | The density of the ice above the basal surface, assumed uniform |
||
| real(kind=r8), | intent(in) | :: | ice_temperature | The temperature of the ice above the basal surface, assumed uniform |
||
| real(kind=r8), | intent(in) | :: | time | The time to which the basal surface should be solved |
||
| logical, | intent(out) | :: | success | True if the solver is successful, false otherwise |
subroutine plume_solve(this, ice_thickness, ice_density, ice_temperature, &
time, success)
!* Author: Chris MacMackin
! Date: March 2017
!
! Solves the state of the plume for the specified ice properties,
! at the specified time. This is done using the a
! quasilinearisation method and a GMRES iterative linear solver.
!
class(plume), intent(inout) :: this
class(scalar_field), intent(in) :: ice_thickness
!! Thickness of the ice above the basal surface
real(r8), intent(in) :: ice_density
!! The density of the ice above the basal surface, assumed uniform
real(r8), intent(in) :: ice_temperature
!! The temperature of the ice above the basal surface, assumed uniform
real(r8), intent(in) :: time
!! The time to which the basal surface should be solved
logical, intent(out) :: success
!! True if the solver is successful, false otherwise
real(r8), dimension(:), allocatable :: solution
real(r8) :: residual
integer, dimension(5) :: info
integer :: flag
class(scalar_field), pointer :: b
call ice_thickness%guard_temp()
b => -ice_thickness/this%r_val
call b%guard_temp()
select type(bound => this%boundaries)
class is(upstream_plume_boundary)
call bound%calculate(time, non_diff_terms, b)
class default
call bound%set_time(time)
end select
solution = this%state_vector()
#ifdef DEBUG
call logger%debug('plume%solve','Calling QLM ODE solver')
#endif
call quasilinear_solve(L, f, jac_prod, solution, 1, residual, flag, info, &
1.e-9_r8*size(solution), precond=preconditioner, &
iter_max=100, krylov_dim=85, gmres_iter_max=5000)
call this%update(solution)
#ifdef DEBUG
call logger%debug('plume%solve','QLM solver required '// &
trim(str(info(5)))//' nonlinear iterations '// &
'and '//trim(str(info(1)+info(2)))// &
' function calls.')
#endif
select case(flag)
case(0)
call logger%trivia('plume%solver','Solved plume at time '//trim(str(time)))
success = .true.
this%time = time
case(1)
call logger%warning('plume%solver','Plume solver stagnated with '// &
'residual of '//trim(str(residual)))
success = .false.
case(2)
call logger%error('plume%solver','Reached maximum number of '// &
'iterations solving plume')
success = .false.
case(3)
call logger%error('plume%solver','Plume solution began to diverge.')
success = .false.
case default
call logger%error('plume%solve','QLM solver failed for plume with '// &
'error code '//trim(str(flag)))
success = .false.
end select
call ice_thickness%clean_temp(); call b%clean_temp()
contains
function L(v)
!! The linear differentiation operator
real(r8), dimension(:), intent(in) :: v
!! The state vector for the system of differential equations
real(r8), dimension(size(v)) :: L
integer :: st, en, btype_l, btype_u, bdepth_l, bdepth_u
type(cheb1d_scalar_field) :: scalar_tmp
type(cheb1d_vector_field) :: vector_tmp
type(cheb1d_vector_field) :: coriolis
call this%update(v)
! Thickness
scalar_tmp = this%thickness%d_dx(1)
call this%boundaries%thickness_bound_info(-1, btype_l, bdepth_l)
call this%boundaries%thickness_bound_info(1, btype_u, bdepth_u)
if (this%lower_bounds(1)) then
call scalar_tmp%set_boundary(-1, 1, this%thickness%get_boundary(-1, 1))
end if
if (this%upper_bounds(1)) then
call scalar_tmp%set_boundary(1, 1, this%thickness%get_boundary(1, 1))
end if
st = 1
en = st + this%thickness_size - 1
L(st:en) = scalar_tmp%raw()
! Velocity
vector_tmp = this%velocity%d_dx(1) - this%velocity_dx
if (this%lower_bounds(2)) then
call vector_tmp%set_boundary(-1, 1, this%velocity%get_boundary(-1, 1))
end if
if (this%upper_bounds(2)) then
call vector_tmp%set_boundary(1, 1, this%velocity%get_boundary(1, 1))
end if
st = en + 1
en = st + this%velocity_size - 1
L(st:en) = vector_tmp%raw()
if (this%phi /= 0._r8) then
coriolis = [0._r8, 0._r8, this%phi/this%nu] .cross. this%velocity
vector_tmp = this%velocity_dx%d_dx(1) - coriolis
else
vector_tmp = this%velocity_dx%d_dx(1)
end if
if (this%lower_bounds(3)) then
call vector_tmp%set_boundary(-1, 1, this%velocity_dx%get_boundary(-1, 1))
end if
if (this%upper_bounds(3)) then
call vector_tmp%set_boundary(1, 1, this%velocity_dx%get_boundary(1, 1))
end if
st = en + 1
en = st + this%velocity_size - 1
L(st:en) = vector_tmp%raw()
! Temperature
scalar_tmp = this%temperature%d_dx(1) - this%temperature_dx
if (this%lower_bounds(4)) then
call scalar_tmp%set_boundary(-1, 1, this%temperature%get_boundary(-1, 1))
end if
if (this%upper_bounds(4)) then
call scalar_tmp%set_boundary(1, 1, this%temperature%get_boundary(1, 1))
end if
st = en + 1
en = st + this%temperature_size - 1
L(st:en) = scalar_tmp%raw()
scalar_tmp = this%temperature_dx%d_dx(1)
if (this%lower_bounds(5)) then
call scalar_tmp%set_boundary(-1, 1, this%temperature_dx%get_boundary(-1, 1))
end if
if (this%upper_bounds(5)) then
call scalar_tmp%set_boundary(1, 1, this%temperature_dx%get_boundary(1, 1))
end if
st = en + 1
en = st + this%temperature_size - 1
L(st:en) = scalar_tmp%raw()
! Salinity
scalar_tmp = this%salinity%d_dx(1) - this%salinity_dx
if (this%lower_bounds(6)) then
call scalar_tmp%set_boundary(-1, 1, this%salinity%get_boundary(-1, 1))
end if
if (this%upper_bounds(6)) then
call scalar_tmp%set_boundary(1, 1, this%salinity%get_boundary(1, 1))
end if
st = en + 1
en = st + this%salinity_size - 1
L(st:en) = scalar_tmp%raw()
scalar_tmp = this%salinity_dx%d_dx(1)
if (this%lower_bounds(7)) then
call scalar_tmp%set_boundary(-1, 1, this%salinity_dx%get_boundary(-1, 1))
end if
if (this%upper_bounds(7)) then
call scalar_tmp%set_boundary(1, 1, this%salinity_dx%get_boundary(1, 1))
end if
st = en + 1
en = st + this%salinity_size - 1
L(st:en) = scalar_tmp%raw()
end function L
subroutine non_diff_terms(D, Uvec, T, S, b, DU_x, DUU_x, DUT_x, DUS_x)
!! Computes the values of \((DU)_x\), \((DU\vec{U})_x\),
!! \((DUT)_x\), \((DUS)_x\), when diffusion is not
!! included. This should be able to handle uniform field types,
!! for use in an ODE solver when integrating near the
!! boundary. The momentum term is calculated quite differently
!! depending on whether it is being done for the ODE solver
!! orthe plume solver. For this reason, as well as the
!! usefulness of avoiding taking any derivatives when
!! calculating it, I have chosen not to bother calculating it
!! here unless it is for the ODE solver.
class(scalar_field), intent(in) :: D
!! The plume thickness
class(vector_field), intent(in) :: Uvec
!! The plume velocity
class(scalar_field), intent(in) :: T
!! The plume temperature
class(scalar_field), intent(in) :: S
!! The plume salinity
class(scalar_field), intent(in) :: b
!! The debth of the base of the ice shelf
class(scalar_field), intent(out) :: DU_x
!! The derivative of the product DU
class(vector_field), intent(out) :: DUU_x
!! The derivative of the product DUU
class(scalar_field), intent(out) :: DUT_x
!! The derivative of the product DUT
class(scalar_field), intent(out) :: DUS_x
!! The derivative of the product DUS
integer :: dims
class(scalar_field), pointer :: m, rho, e, S_a, U, V, &
T_a, rho_a, rho_x, Unorm
class(scalar_field), allocatable, dimension(:) :: tmp
type(cheb1d_vector_field) :: coriolis
call D%guard_temp(); call Uvec%guard_temp(); call T%guard_temp()
call S%guard_temp(); call b%guard_temp()
S_a => this%ambient_conds%ambient_salinity(b,time)
T_a => this%ambient_conds%ambient_temperature(b,time)
call S_a%guard_temp(); call T_a%guard_temp()
rho => this%eos%water_density(T, S)
rho_a => this%eos%water_density(T_a, S_a)
U => Uvec%component(1)
V => Uvec%component(2)
call rho%guard_temp(); call rho_a%guard_temp(); call U%guard_temp()
call V%guard_temp()
e => this%entrainment_formulation%entrainment_rate(Uvec, D, b, rho_a - rho, time)
call e%guard_temp()
call this%melt_formulation%solve_for_melt(Uvec, b, T, S, D, time)
m => this%melt_formulation%melt_rate()
call m%guard_temp()
DU_x = e + m
if (this%melt_formulation%has_heat_terms()) then
DUT_x = e*T_a - this%melt_formulation%heat_equation_terms()
else
DUT_x = e*T_a
end if
if (this%melt_formulation%has_salt_terms()) then
DUS_x = e*S_a - this%melt_formulation%salt_equation_terms()
! print*, DUS_x%raw()
else
DUS_x = e*S_a
end if
select type(Uvec)
class is(uniform_vector_field)
Unorm => Uvec%norm()
rho_x => this%eos%water_density_derivative(T, (DUT_x - DU_x*T)/(D*U), &
S, (DUS_x - DU_x*S)/(D*U), 1)
call Unorm%guard_temp(); call rho_x%guard_temp()
dims = Uvec%raw_size()/Uvec%elements()
allocate(tmp(dims), mold=D)
tmp(1) = 1._r8 - this%delta*D*(rho_a - rho)/U**2
!print*,tmp(1)%raw()
tmp(1) = (D*(rho_a - rho)*(b%d_dx(1) - 2*this%delta*DU_x/U) &
+ 0.5*this%delta*D**2*rho_x - this%mu*Unorm*U &
+ this%phi*D*V) / (1._r8 - this%delta*D*(rho_a - rho)/U**2)
if (dims > 1) then
tmp(2) = -this%mu*Unorm*V - this%phi*D*U
end if
DUU_x = tmp
call Unorm%clean_temp(); call rho_x%clean_temp()
class default
! if (this%phi /= 0._r8) then
! coriolis = [0._r8, 0._r8, this%phi] .cross. Uvec
! DUU_x = -this%mu*Uvec*Uvec%norm() + 0.5_r8*this%delta*D**2*(.grad. rho) &
! + D*(rho_a - rho)*(.grad.(b - this%delta*D)) - D*coriolis
! else
DUU_x = -this%mu*Uvec*Uvec%norm() + 0.5_r8*this%delta*D**2*(.grad. rho) &
+ D*(rho_a - rho)*(.grad.(b - this%delta*D))
! end if
end select
call e%clean_temp(); call S_a%clean_temp(); call T_a%clean_temp()
call rho%clean_temp(); call m%clean_temp(); call rho_a%clean_temp()
call U%clean_temp(); call V%clean_temp()
call D%clean_temp(); call Uvec%clean_temp(); call T%clean_temp()
call S%clean_temp(); call b%clean_temp()
end subroutine non_diff_terms
function f(v)
!! The nonlinear operator
real(r8), dimension(:,:), intent(in) :: v
!! The state vector for the system of differential equations,
!! and its derivatives. Column \(i\) represents the \(i-1\)
!! derivative.
real(r8), dimension(size(v,1)) :: f
call this%update(v(:,1))
call nonlinear(f, .false.)
end function f
function jac_prod(v, dv)
!! The product of the Jacobian of the nonlienar operator at v,
!! multiplying dv.
real(r8), dimension(:,:), intent(in) :: v
!! The state vector for the system of differential equations,
!! and its derivatives. Column \(i\) represents the \(i-1\)
!! derivative.
real(r8), dimension(:,:), intent(in) :: dv
!! The state vector for the system of differential equations,
!! and its derivatives, to be multiplied by the
!! Jacobian. Column \(i\) represents the \(i-1\) derivative.
real(r8), dimension(size(v,1)) :: jac_prod
type(cheb1d_scalar_field) :: stmp
type(cheb1d_vector_field) :: vtmp
integer :: i
call this%update(v(:,1))
call stmp%assign_meta_data(this%thickness)
call vtmp%assign_meta_data(this%velocity)
call stmp%set_from_raw(dv(1:this%thickness_size, 1))
call this%thickness%set_derivative(stmp)
i = 1 + this%thickness_size
call vtmp%set_from_raw(dv(i:i + this%velocity_size - 1, 1))
call this%velocity%set_derivative(vtmp)
i = i + this%velocity_size
call vtmp%set_from_raw(dv(i:i + this%velocity_size - 1, 1))
call this%velocity_dx%set_derivative(vtmp)
i = i + this%velocity_size
call stmp%set_from_raw(dv(i:i + this%temperature_size - 1, 1))
call this%temperature%set_derivative(stmp)
i = i + this%temperature_size
call stmp%set_from_raw(dv(i:i + this%temperature_size - 1, 1))
call this%temperature_dx%set_derivative(stmp)
i = i + this%temperature_size
call stmp%set_from_raw(dv(i:i + this%salinity_size - 1, 1))
call this%salinity%set_derivative(stmp)
i = i + this%salinity_size
call stmp%set_from_raw(dv(i:i + this%salinity_size - 1, 1))
call this%salinity_dx%set_derivative(stmp)
call nonlinear(jac_prod, .true.)
call this%thickness%unset_derivative()
call this%velocity%unset_derivative()
call this%velocity_dx%unset_derivative()
call this%temperature%unset_derivative()
call this%temperature_dx%unset_derivative()
call this%salinity%unset_derivative()
call this%salinity_dx%unset_derivative()
end function jac_prod
subroutine nonlinear(f, deriv)
real(r8), dimension(:), intent(out) :: f
logical, intent(in) :: deriv
!! If true, return Jacobian product, otherwise return result
!! of nonlienar operator.
integer :: st, en
type(cheb1d_scalar_field) :: scalar_tmp, D_x, D_nd, S_nd, T_nd
type(cheb1d_vector_field) :: vector_tmp, U_nd
class(scalar_field), pointer :: U, U_x
! Use same or similar notation for variables as used in equations
associate(D => this%thickness, Uvec => this%velocity, &
Uvec_x => this%velocity_dx, S => this%salinity, &
S_x => this%salinity_dx, T => this%temperature, &
T_x => this%temperature_dx, mf => this%melt_formulation, &
h => ice_thickness, delta => this%delta, nu => this%nu, &
mu => this%mu, r => this%r_val, bounds => this%boundaries)
call non_diff_terms(D, Uvec, T, S, b, D_nd, U_nd, T_nd, S_nd)
! FIXME: Alter this so that can take advantage of
! parameterisations returning uniform fields.
U => this%velocity%component(1)
U_x => this%velocity_dx%component(1)
call U%guard_temp(); call U_x%guard_temp()
! Thickness
scalar_tmp = (D_nd - D*U_x)/U
D_x = scalar_tmp
if (this%lower_bounds(1)) then
call scalar_tmp%set_boundary(-1, 1, bounds%thickness_bound(-1))
end if
if (this%upper_bounds(1)) then
call scalar_tmp%set_boundary(1, 1, bounds%thickness_bound(1))
end if
st = 1
en = st + this%thickness_size - 1
if (deriv) then
scalar_tmp = scalar_tmp%get_derivative()
end if
f(st:en) = scalar_tmp%raw()
! Velocity
vector_tmp = 0._r8 * Uvec
if (this%lower_bounds(2)) then
call vector_tmp%set_boundary(-1, 1, bounds%velocity_bound(-1))
end if
if (this%upper_bounds(2)) then
call vector_tmp%set_boundary(1, 1, bounds%velocity_bound(1))
end if
st = en + 1
en = st + this%velocity_size - 1
if (deriv) then
vector_tmp = vector_tmp%get_derivative()
end if
f(st:en) = vector_tmp%raw()
vector_tmp = D*U*Uvec_x !Needed due to compiler bug
vector_tmp = (vector_tmp + D*U_x*Uvec + D_x*U*Uvec - U_nd - &
nu*D_x*Uvec_x)/(nu*D)
if (this%lower_bounds(3)) then
call vector_tmp%set_boundary(-1, 1, bounds%velocity_bound(-1))
end if
if (this%upper_bounds(3)) then
call vector_tmp%set_boundary(1, 1, bounds%velocity_bound(1))
end if
st = en + 1
en = st + this%velocity_size - 1
if (deriv) then
vector_tmp = vector_tmp%get_derivative()
end if
f(st:en) = vector_tmp%raw()
! Temperature
scalar_tmp = uniform_scalar_field(0._r8)
if (this%lower_bounds(4)) then
call scalar_tmp%set_boundary(-1, 1, bounds%temperature_bound(-1))
end if
if (this%upper_bounds(4)) then
call scalar_tmp%set_boundary(1, 1, bounds%temperature_bound(1))
end if
st = en + 1
en = st + this%temperature_size - 1
if (deriv) then
scalar_tmp = scalar_tmp%get_derivative()
end if
f(st:en) = scalar_tmp%raw()
scalar_tmp = (D*U*T_x + D*U_x*T + D_x*U*T - T_nd - nu*D_x*T_x)/(nu*D)
if (this%lower_bounds(5)) then
call scalar_tmp%set_boundary(-1, 1, bounds%temperature_bound(-1))
end if
if (this%upper_bounds(5)) then
call scalar_tmp%set_boundary(1, 1, bounds%temperature_bound(1))
end if
st = en + 1
en = st + this%salinity_size - 1
if (deriv) then
scalar_tmp = scalar_tmp%get_derivative()
end if
f(st:en) = scalar_tmp%raw()
! Salinity
scalar_tmp = uniform_scalar_field(0._r8)
if (this%lower_bounds(6)) then
call scalar_tmp%set_boundary(-1, 1, bounds%salinity_bound(-1))
end if
if (this%upper_bounds(6)) then
call scalar_tmp%set_boundary(1, 1, bounds%salinity_bound(1))
end if
st = en + 1
en = st + this%salinity_size - 1
if (deriv) then
scalar_tmp = scalar_tmp%get_derivative()
end if
f(st:en) = scalar_tmp%raw()
scalar_tmp = (D*U*S_x + D*U_x*S + D_x*U*S - S_nd - nu*D_x*S_x)/(nu*D)
if (this%lower_bounds(7)) then
call scalar_tmp%set_boundary(-1, 1, bounds%salinity_bound(-1))
end if
if (this%upper_bounds(7)) then
call scalar_tmp%set_boundary(1, 1, bounds%salinity_bound(1))
end if
st = en + 1
en = st + this%salinity_size - 1
if (deriv) then
scalar_tmp = scalar_tmp%get_derivative()
end if
f(st:en) = scalar_tmp%raw()
call U%clean_temp(); call U_x%clean_temp()
end associate
end subroutine nonlinear
function preconditioner(v, state, L_op, f_op, fcur, rhs)
!! The preconditioner, which approximates an inverse of `L`.
real(r8), dimension(:), intent(in) :: v
!! The vector to be preconditioned.
real(r8), dimension(:,:), intent(in) :: state
!! The current state vector for the system of differential
!! equations, and its derivatives. Column \(i\) represents the
!! \(i-1\) derivative.
procedure(L) :: L_op
!! The linear, left-hand-side of the ODE being solved.
procedure(f) :: f_op
!! The nonlinear, right-hand-side of the ODE being solved.
real(r8), dimension(:), intent(in) :: fcur
!! The result of `f(u)`
real(r8), dimension(:), intent(in) :: rhs
!! The right hand side of the linear system being
!! preconditioned.
real(r8), dimension(size(v)) :: preconditioner
!! The result of applying the preconditioner.
integer :: st, en, ust, uen, pst, pen
integer :: bloc
real(r8) :: nu
type(plume) :: v_plume
type(cheb1d_scalar_field) :: scalar_tmp
type(cheb1d_vector_field) :: vector_tmp
class(scalar_field), pointer :: U, U_x
v_plume%thickness_size = this%thickness_size
call v_plume%thickness%assign_meta_data(this%thickness)
v_plume%velocity_size = this%velocity_size
call v_plume%velocity%assign_meta_data(this%velocity)
call v_plume%velocity_dx%assign_meta_data(this%velocity_dx)
v_plume%temperature_size = this%temperature_size
call v_plume%temperature%assign_meta_data(this%temperature)
call v_plume%temperature_dx%assign_meta_data(this%temperature_dx)
v_plume%salinity_size = this%salinity_size
call v_plume%salinity%assign_meta_data(this%salinity)
call v_plume%salinity_dx%assign_meta_data(this%salinity_dx)
call v_plume%update(v)
nu = this%nu
bloc = get_bound_loc(1)
v_plume%thickness = this%precond%solve_for(v_plume%thickness, bloc, &
v_plume%thickness%get_boundary(bloc, 1))
st = 1
en = st + this%thickness_size - 1
preconditioner(st:en) = v_plume%thickness%raw()
! Precondition the U_x term after have preconditioned values for S and T
st = en + 1
en = st + this%velocity_size - 1
ust = st
uen = en
bloc = get_bound_loc(3)
v_plume%velocity_dx = this%precond%solve_for(v_plume%velocity_dx, bloc, &
v_plume%velocity_dx%get_boundary(bloc, 1))
st = en + 1
en = st + this%velocity_size - 1
preconditioner(st:en) = v_plume%velocity_dx%raw()
bloc = get_bound_loc(2)
vector_tmp = v_plume%velocity + v_plume%velocity_dx
v_plume%velocity = this%precond%solve_for(vector_tmp, bloc, &
v_plume%velocity%get_boundary(bloc, 1))
preconditioner(ust:uen) = v_plume%velocity%raw()
! Precondition T_x terms before T
st = en + 1
en = st + this%temperature_size - 1
pst = st
pen = en
bloc = get_bound_loc(5)
v_plume%temperature_dx = this%precond%solve_for(v_plume%temperature_dx, bloc, &
v_plume%temperature_dx%get_boundary(bloc, 1))
st = en + 1
en = st + this%temperature_size - 1
preconditioner(st:en) = v_plume%temperature_dx%raw()
bloc = get_bound_loc(4)
v_plume%temperature = this%precond%solve_for(v_plume%temperature + &
v_plume%temperature_dx, bloc, v_plume%temperature%get_boundary(bloc, 1))
preconditioner(pst:pen) = v_plume%temperature%raw()
! Precondition S_x terms before S
st = en + 1
en = st + this%salinity_size - 1
pst = st
pen = en
bloc = get_bound_loc(7)
v_plume%salinity_dx = this%precond%solve_for(v_plume%salinity_dx, bloc, &
v_plume%salinity_dx%get_boundary(bloc, 1))
st = en + 1
en = st + this%temperature_size - 1
preconditioner(st:en) = v_plume%salinity_dx%raw()
bloc = get_bound_loc(6)
v_plume%salinity = this%precond%solve_for(v_plume%salinity+v_plume%salinity_dx, bloc, &
v_plume%salinity%get_boundary(bloc, 1))
preconditioner(pst:pen) = v_plume%salinity%raw()
end function preconditioner
integer function get_bound_loc(component_id)
integer :: component_id
if (this%lower_bounds(component_id)) then
get_bound_loc = -1
else if (this%upper_bounds(component_id)) then
get_bound_loc = 1
else
get_bound_loc = 0
end if
end function get_bound_loc
end subroutine plume_solve