I'm trying to use the output temperature from HeatTransportSolver to update the fluidity field passed to FlowSolver defined on a hybrid model with the same geometry. I get an error when passing the result of icepack.rate_factor(T). This doesn't happen when I calculate the rate factor manually, i.e using A = firedrake.interpolate(firedrake.exp(-60.0/(ideal_gas*T)),Q)
import icepack
import firedrake
from icepack.constants import (
year,
thermal_diffusivity as alpha,
ideal_gas
)
"""
Set up same problem as in heat_transport_test.py
"""
Nx, Ny = 32, 1
Lx, Ly = 20e3, 1e3
mesh2d = firedrake.RectangleMesh(Nx, Ny, Lx, Ly)
mesh = firedrake.ExtrudedMesh(mesh2d, layers=1)
Q = firedrake.FunctionSpace(mesh, "CG", 2, vfamily="GL", vdegree=4)
Q_c = firedrake.FunctionSpace(mesh, "CG", 2, vfamily="R", vdegree=0)
V = firedrake.VectorFunctionSpace(mesh, "CG", 2, vfamily="GL", vdegree=4, dim=2)
W = firedrake.FunctionSpace(mesh, "DG", 1, vfamily="GL", vdegree=5)
x, y, z = firedrake.SpatialCoordinate(mesh)
h0, dh = 500.0, 100.0
h = firedrake.interpolate(h0 - dh * x / Lx, Q_c)
s0, ds = 500.0, 50.0
s = firedrake.interpolate(s0 - ds * x / Lx, Q_c)
E_surface = 480
q_bed = 50e-3 * year * 1e-6
E_initial = firedrake.interpolate(E_surface + q_bed / alpha * h * (1 - z), Q)
u0 = 100.0
du = 100.0
u_expr = firedrake.as_vector((u0 + du * x / Lx, 0))
u = firedrake.interpolate(u_expr, V)
w = firedrake.interpolate((-du / Lx + dh / Lx / h * u[0]) * z, W)
E_q = E_initial.copy(deepcopy=True)
E_0 = E_initial.copy(deepcopy=True)
heat_model = icepack.models.HeatTransport3D()
T = heat_model.temperature(E_q)
A = icepack.rate_factor(T)
fields = {
"velocity": u,
"vertical_velocity": w,
"thickness": h,
"surface": s,
"heat_bed": firedrake.Constant(q_bed),
"energy_inflow": E_initial,
"energy_surface": firedrake.Constant(E_surface),
}
dt = 10.0
final_time = Lx / u0
num_steps = int(final_time / dt) + 1
heat_solver = icepack.solvers.HeatTransportSolver(heat_model,dirichlet_ids=[1],side_wall_ids = [3,4])
"""
Now define a Hybrid model over the same domain to solve for the flow.
"""
hybrid_model = icepack.models.HybridModel()
flow_solver = icepack.solvers.FlowSolver(hybrid_model, dirichlet_ids=[1],side_wall_ids=[3,4])
#Solve for energy
E = heat_solver.solve(
dt,
energy=E_0,
velocity=u,
vertical_velocity=w,
thickness=h,
surface=s,
heat=firedrake.Constant(0.0),
heat_bed=firedrake.Constant(q_bed),
energy_inflow=E_initial,
energy_surface=firedrake.Constant(E_surface),
)
#obtain temperature field
T = heat_model.temperature(E)
#calculate rate factor field
A = firedrake.interpolate(firedrake.exp(-60.0/(ideal_gas*T)),Q) #this works!
#A = icepack.rate_factor(T) #this doesn't
#run flow solver using rate factor field
u0 = flow_solver.diagnostic_solve(
velocity=u,
thickness=h,
surface=s,
#fluidity=icepack.rate_factor(255.0),
fluidity =A,
friction = firedrake.Constant(0.0)
)
TypeError Traceback (most recent call last)
/var/folders/4q/8fx0bbzx37x2pnkf3h36dsgm0000gn/T/ipykernel_5277/1250458578.py in <module>
93
94 #run flow solver using rate factor field
---> 95 u0 = flow_solver.diagnostic_solve(
96 velocity=u,
97 thickness=h,
~/Documents/Python_Envs/icepack/icepack/solvers/flow_solver.py in diagnostic_solve(self, **kwargs)
136 def diagnostic_solve(self, **kwargs):
137 r"""Solve the diagnostic model physics for the ice velocity"""
--> 138 return self._diagnostic_solver.solve(**kwargs)
139
140 def prognostic_solve(self, dt, **kwargs):
~/Documents/Python_Envs/icepack/icepack/solvers/flow_solver.py in solve(self, **kwargs)
208 r"""Solve the diagnostic model physics for the ice velocity"""
209 if not hasattr(self, "_solver"):
--> 210 self.setup(**kwargs)
211 else:
212 for name, field in kwargs.items():
~/Documents/Python_Envs/icepack/icepack/solvers/flow_solver.py in setup(self, **kwargs)
168 self._fields[name] = field.copy(deepcopy=True)
169 else:
--> 170 raise TypeError(
171 "Input %s field has type %s, must be Constant or Function!"
172 % (name, type(field))
TypeError: Input fluidity field has type <class 'ufl.algebra.Product'>, must be Constant or Function!