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(asym_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 asym_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(asym_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('asym_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('asym_plume%solver','Solved plume at time '//trim(str(time)))
success = .true.
this%time = time
case(1)
call logger%warning('asym_plume%solver','Plume solver stagnated with '// &
'residual of '//trim(str(residual)))
success = .false.
case(2)
call logger%error('asym_plume%solver','Reached maximum number of '// &
'iterations solving plume')
success = .false.
case(3)
call logger%error('asym_plume%solver','Plume solution began to diverge.')
success = .false.
case default
call logger%error('asym_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, ctmp(2)
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
ctmp(1) = -this%phi*this%shape%a_DV/(this%nu*this%shape%a_DU) &
*this%velocity%component(2)
ctmp(2) = this%phi*this%shape%a_DU/(this%nu*this%shape%a_DV) &
*this%velocity%component(1)
coriolis = ctmp
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()
! print*,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.
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
class(scalar_field), pointer :: m, rho_b, rho_t, 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)
select type(eos => this%eos)
class is(ave_linear_eos)
rho_b => eos%water_density_ave1(T, S)
rho_t => eos%water_density_ave2(T, S)
class default
rho_b => eos%water_density(T, S)
rho_t => rho_b
end select
U => Uvec%component(1)
V => Uvec%component(2)
call S_a%guard_temp(); call T_a%guard_temp(); call rho_b%guard_temp()
call rho_t%guard_temp(); call U%guard_temp(); call V%guard_temp()
rho_a => this%eos%water_density(T_a, S_a)
call rho_a%guard_temp()
e => this%entrainment_formulation%entrainment_rate(Uvec, D, b, rho_a - rho_b, 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()
associate(f_D1 => this%shape%f_D1, f_D2 => this%shape%f_D2, &
f_U2 => this%shape%f_U2, f_V2 => this%shape%f_V2, &
f_S1 => this%shape%f_S1, f_S2 => this%shape%f_S2, &
f_T1 => this%shape%f_T1, f_T2 => this%shape%f_T2, &
a_DU => this%shape%a_DU, a_DU2 => this%shape%a_DU2, &
a_DUV => this%shape%a_DUV, a_DUS => this%shape%a_DUS, &
a_DUT => this%shape%a_DUT, a_D2 => this%shape%a_D2, &
a_DV => this%shape%a_DV, a_UabsU => this%shape%a_UabsU, &
a_UabsV => this%shape%a_UabsV, dy => this%dy)
DU_x = (e + m - f_D2*D * f_V2*V/dy)/a_DU
if (this%melt_formulation%has_heat_terms()) then
DUT_x = (e*T_a - f_D2*D * f_V2*V * f_T2*T/dy - &
this%melt_formulation%heat_equation_terms())/a_DUT
else
DUT_x = (e*T_a - f_D2*D * f_V2*V * f_T2*T/dy)/a_DUT
end if
if (this%melt_formulation%has_salt_terms()) then
DUS_x = (e*S_a - f_D2*D * f_V2*V * f_S2*S/dy - &
this%melt_formulation%salt_equation_terms())/a_DUS
else
DUS_x = (e*S_a - f_D2*D * f_V2*V * f_S2*S/dy)/a_DUS
end if
Unorm => Uvec%norm()
call Unorm%guard_temp()
select type(Uvec)
class is(uniform_vector_field)
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 rho_x%guard_temp()
allocate(tmp(this%vel_dims), mold=D)
tmp(1) = (D*(rho_a - rho_b)*b%d_dx(1) &
- 2._r8*D*this%delta*a_D2*(rho_a - rho_t)*DU_x/U &
+ 0.5*this%delta*a_D2*D**2*rho_x &
- this%mu*a_UabsU*Unorm*U &
+ this%phi*a_DV*D*V &
- f_D2*D * f_U2*U * f_V2*V/dy) / &
(a_DU2 - this%delta*a_D2*D*(rho_a - rho_t)/U**2)
if (this%vel_dims > 1) then
! Use entrainment and melt as work-arrays to hold the
! upper and lower boundary density values. Should be able
! to just multiply the function results directly, but
! there's a compiler bug.
call e%clean_temp(); call m%clean_temp()
e => this%eos%water_density(f_T2*T, f_S2*S)
m => this%eos%water_density(f_T1*T, f_S1*S)
call e%guard_temp(); call m%guard_temp()
tmp(2) = ((f_D2**2*(e - rho_a) - f_D1**2*(m - rho_a))* &
0.5_r8*this%delta/dy*D**2 - this%mu*a_UabsV*Unorm*V &
- a_DU*this%phi*D*U - f_D2*D * f_V2**2*V**2/dy)/a_DUV
end if
DUU_x = tmp
call rho_x%clean_temp()
class default
allocate(tmp(this%vel_dims), mold=D)
tmp(1) = b%d_dx(1)
tmp(1) = (D*(rho_a - rho_b)*b%d_dx(1) &
- a_D2*D*(rho_a - rho_t)*this%delta*D%d_dx(1) &
+ 0.5_r8*a_D2*this%delta*D**2*rho_t%d_dx(1) &
- this%mu*a_UabsU*U*Unorm &
- f_D2*D * f_V2*V * f_U2*U/dy)/a_DU2
if (this%vel_dims > 1) then
! Use entrainment and melt as work-arrays to hold the
! upper and lower boundary density values. Should be able
! to just multiply the function results directly, but
! there's a compiler bug.
call e%clean_temp(); call m%clean_temp()
e => this%eos%water_density(f_T2*T, f_S2*S)
m => this%eos%water_density(f_T1*T, f_S1*S)
call e%guard_temp(); call m%guard_temp()
tmp(2) = ((f_D2**2*(e - rho_a) - f_D1**2*(m - rho_a))* &
0.5_r8*this%delta/dy*D**2 - this%mu*a_UabsV*V*Unorm &
- f_D2*D * f_V2**2 * V**2/dy)/a_DUV
end if
DUU_x = tmp
end select
call Unorm%clean_temp()
end associate
call e%clean_temp(); call S_a%clean_temp(); call T_a%clean_temp()
call rho_b%clean_temp(); call m%clean_temp(); call rho_a%clean_temp()
call rho_t%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_scalar_field), allocatable, dimension(:) :: vtmp
type(cheb1d_vector_field) :: vector_tmp, U_nd
class(scalar_field), pointer :: U, V, U_x, V_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,&
a_DU2 => this%shape%a_DU2, a_DUV => this%shape%a_DUV, &
a_DUS => this%shape%a_DUS, a_DUT => this%shape%a_DUT, &
a_DU => this%shape%a_DU, a_DV => this%shape%a_DV, &
a_DS => this%shape%a_DS, a_DT => this%shape%a_DT, &
f_D2 => this%shape%f_D2, f_Up => this%shape%f_Up, &
f_Vp => this%shape%f_Vp, f_Sp => this%shape%f_Sp, &
f_Tp => this%shape%f_Tp, dy => this%dy )
call non_diff_terms(D, Uvec, T, S, b, D_nd, U_nd, T_nd, S_nd)
U => this%velocity%component(1)
U_x => this%velocity_dx%component(1)
call U%guard_temp(); call U_x%guard_temp()
if (this%vel_dims > 1) then
V => this%velocity%component(2)
V_x => this%velocity_dx%component(2)
call V%guard_temp(); call V_x%guard_temp()
end if
! 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()
allocate (vtmp(this%vel_dims), mold=D)
vtmp(1) = (a_DU2*(2._r8*D*U*U_x + D_x*U**2 - U_nd%component(1)) &
- nu*a_DU*D_x*U_x - nu*f_D2*D*f_Up*U/dy)/(nu*a_DU*D)
if (this%vel_dims > 1) then
vtmp(2) = (a_DUV*(D*U*V_x + D*U_x*V + D_x*U*V - U_nd%component(2)) &
- nu*a_DV*D_x*V_x - nu*f_D2*D*f_Vp*V/dy)/(nu*a_DV*D)
end if
vector_tmp = vtmp
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 = (a_DUT*(D*U*T_x + D*U_x*T + D_x*U*T - T_nd) - &
nu*a_DT*D_x*T_x - nu*f_D2*D*f_Tp*T/dy)/(nu*a_DT*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 = (a_DUS*(D*U*S_x + D*U_x*S + D_x*U*S - S_nd) - &
nu*a_DS*D_x*S_x - nu*f_D2*D*f_Sp*S/dy)/(nu*a_DS*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()
if (this%vel_dims > 0) then
call V%clean_temp(); call V_x%clean_temp()
end if
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, i
real(r8) :: nu
type(asym_plume) :: v_plume
type(cheb1d_scalar_field) :: scalar_tmp
type(cheb1d_vector_field) :: vector_tmp, tmp2
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(1, 1), -1)
st = 1
en = st + this%thickness_size - 1
preconditioner(st:en) = v_plume%thickness%raw()
bloc = get_bound_loc(3)
v_plume%velocity_dx = this%precond%solve_for(v_plume%velocity_dx, bloc, &
v_plume%velocity_dx%get_boundary(1, 1), -1)
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(1, 1), -1)
st = en + 1
en = st + this%velocity_size - 1
preconditioner(st:en) = v_plume%velocity%raw()
st = en + 1
en = st + this%velocity_size - 1
preconditioner(st:en) = v_plume%velocity_dx%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(1, 1), -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(1, 1), -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(1, 1), -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(1, 1), -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 asym_plume_solve