icepack/icepack

Finite element modeling of glaciers and ice sheets

icepack's Issues

passing expressions to solver objects

The FlowSolver class contains a member called fields, a dictionary which stores the names and values of all the input parameters to the diagnostic and prognostic solver routines. When you pass in a (keyword) argument to either the diagnostic or prognostic solve routines, it gets assigned into the corresponding copy in the fields dictionary. If the argument is a Function or a Constant, everything is fine.

The problem happens if you pass in a symbolic expression, something more complicated than a Function or Constant like A_0 * exp(-Q / (R * T)) where T is a Function. Expressions can't be assigned to -- changing an expression means the C code under the hood has to get recompiled. Right now, the FlowSolver class just skips over these silently. This leads to the very confusing behavior that if you pass a different value on a repeated solve, nothing gets changed and the old value will be reused.

There are a few possible solutions, in order from most to least conservative:

1. Allow it and throw a warning.
2. Figure out what the value shape of that expression is and project it into a suitable function space. Then throw a warning.
3. Disallow passing non-trivial expressions -- you can only pass a Function or Constant. If you want to pass a more complicated expression, either project it to a function space first, or modify the physics to include that expression and pass its coefficients by name as Functions or Constants.

If the icepack paper were out and there were a bunch of users, option 1 would be the only way to go. At this juncture, we could get away with option 3 and possibly save some aggravation down the line.

Negative rate factors with reasonable temperatures

I am having issues with negative rate factors, and thus negative fluidities and viscosities, when working with somewhat complicated temperature profiles. I can try to work out a test case if needed, but this is the gist of what I get:

Perhaps the question is just whether it is safe to do something like A = max_value(small_number, A) in the rate_factor function or whether a better fix is needed.

I think this is happening as a result of a curving temperature profile (i.e. the coldest bit is at depth, and then when the rate factor gets projected onto the function space I somehow get negative values for the fluidity at the surface. The problem is reduced by increasing the vertical order of the function space for the rate factor, but it is not eliminated even at order 8.

Unexpected differences with Hybrid model compared to Elmer/Ice with Cuffey and Paterson A

I am hoping I am just missing something silly here, but I am getting a substantial difference in velocity when using the Icepack hybrid model versus Elmer/Ice on test cases where I think there should be little or no difference. The SIA model in icepack matches Elmer.

The difference is caused by the rate factor--dividing the default value in icepack by a factor of ~2 ** 2.5 seems to get things almost identical. I have attached a minimal working example where the difference is seen; I am simply using built in rate_factor in icepack and I am using the same values in Elmer/Ice (or, rather, their suggested values but I have separately confirmed that those match Cuffey and Paterson). A reasonably recent Elmer/Ice is probably required for the viscosity formulation, but not bleeding edge and no need for any add-ins other than the elmerice.

For some reason I have been having problems with point evaluation in firedrake recently--for that reason, the last two plots in the notebook are not as clear as the could be, but I think it shows things pretty clearly anyway. Hopefully @danshapero can help point out if I am doing something wrong or if this is a bug in the hybrid solver.
mwe.tar.gz

Zero regularization in 05-ice-shelf-inverse tutorial

In the tutorial notebook 05-ice-shelf-inverse.ipynb the variable Θ is redefined to be a Constant, and the return value is zero.

Requirements not installed

On clean installs, requirements.txt is ignored if icepack is installed with firedrake. I don't think this happens with pip/from source. But if a new user follows instructions on https://icepack.github.io/install/ the test will fail and they will be missing:

geojson, gmsh, meshio>=3.3.1, MeshPy, netCDF4, pooch>=1.0.0, pygmsh<=6.1.1, scipy, shapely, tqdm, and rasterio.

I did not fix this, since I am not sure if there is a better solution than just adding
pip install -r requirements.txt
before the pytest line on the instructions.

Evaluating velocities at depth with the hybrid model

I would like to evaluate a velocity field defined in xyz/xz space at an arbitrary array of points within the modeling domain. However, I am a bit confused with the firedrake conventions when evaluating function spaces in the z-coordinate.

For example, in the hybrid flowband tutorial here, there is a description for resolving shear flow using a weighted depth average. Since the use of Lagrange polynomials as weighting functions was confusing to me, I naively thought I could just evaluate the velocity as:

Nx = 500
Nz = 100

# Create list of points in xz-space to evaluate velocity field
x_grid = np.linspace(5.0e3, 30.0e3, Nx)
z_grid = np.linspace(0.0, 1.0, Nz)
Xg, Zg = np.meshgrid(x_grid, z_grid)
points = np.column_stack((Xg.ravel(), Zg.ravel()))
points = [tuple(pt) for pt in points]

# Run evaluation
u_grid = np.array(u(points)).reshape(Xg.shape)

# Make plot
fig, axes = plt.subplots()
axes.set_title("Shear component of ice velocity")
axes.set_ylabel("velocity (m / yr)")
firedrake.plot(u_shear, axes=axes);        # computed earlier in the notebook
axes.plot(x_grid, u_grid[-1] - u_grid[0])  # surface velocity minus base velocity

However, I find the evaluated velocities give me a (surface - base) difference that's about 3x larger than u_shear. Also, I found that firedrake allows me to evaluate z_grid in the domain [-0.4, 1.4], which was surprising as I figured that $\zeta$ was only defined in [0, 1].

Is there something I'm missing here? Is there a better way to reconstruct the velocities at depth?

No changes when running diagnostic_solve

Hello!

When I ran diagnostic solve in ice stream solver whose model is the Hybrid model with an initial velocity from observing dataset, nothing happened, and the solving result was the same as initial velocity data after simulating several timesteps. As a result, abnormal ice thickness occurs in some areas in the shape.

The code for building the model is as follows:

def friction(**kwargs):
    u = kwargs['velocity']
    h = kwargs['thickness']
    s = kwargs['surface']
    C = kwargs['friction']

    return icepack.models.hybrid.bed_friction(
        velocity=u,
        friction=C
    )

model = icepack.models.HybridModel(friction=friction)
opts = {'dirichlet_ids': boundary_ids}
solver = icepack.solvers.FlowSolver(model, **opts)

u = solver.diagnostic_solve(
        velocity=u0,
        thickness=h,
        surface=s,
        fluidity=A,
        friction=C
    )

δt = 1

for t in trange(85):
    h = solver.prognostic_solve(
            δt,
            thickness=h,
            velocity=u,
            accumulation=a,
            # thickness_inflow=h
        )

    s = icepack.compute_surface(thickness=h, bed=bed)
    
    # where the A and C update with timestep as well, the calculation process is omitted here
    u = solver.diagnostic_solve(
            velocity=u0,
            thickness=h,
            surface=s,
            fluidity=A,
            friction=C
        )

    h.interpolate(max_value(h, 0))
    s = icepack.compute_surface(thickness=h, bed=bed)

Is there any issue? Thank you.

separate models and solvers

Right now the model classes use firedrake.solve instead of constructing reusable variational problem / solver objects. This is wasteful computationally, especially for the hybrid model.

On top of that, the model class is responsible both for specifying what problem you're solving as well as how you're solving it. By separating out the "how" into a separate class we can both eliminate the waste and improve the overall design.

The usage will look something like this:

model = icepack.models.IceShelf()
solver = model.create_solver(u=u, h=h, A=A, **opts)
for step in range(num_timesteps):
    u = solver.diagnostic_solve(u=u, h=h, A=A)
    h = solver.prognostic_solve(dt, h=h, u=u, a=a, h_inflow=h_inflow)

There's one extra step, creating the solver object. There is one possible downside though. With the current design, you can do crazy things like using a P1 discretization for the thickness at one step and a P2 discretization at the next step. If all the input fields are getting cached, this would probably be an error. That's pretty obvious to the user, but there could be subtler issues of reference to out-of-date objects silently giving wrong results.

@hoffmaao , @jabadge thoughts?

Error with multi-parameter optimization

icepack/src/icepack/statistics.py

Line 164 in 0934317

self._controls = [field.copy(deepcopy=True) for field in problem.controls]

problem.controls should be self.problem.controls. After this change multi-parameter Adjoint optimization works as expected.

Questions about the exact equation in the Synthetic ice stream (xy) case

Hi, may I ask what will be the exact equation to generate the data in the tutorial of Synthetic ice stream (xy) case?

more specific errors

It would be great to make some of the icepack errors more specific. For example, when I get the error: "TypeError: Input fields must be Constant or Function!", it would be nice to know which fields are not meeting this criteria.

Parallelism in icepack and Firedrake

What is the canonical way to exploit parallelism in icepack/Firedrake/PETSc? Both the Firedrake and PETSc documentation recommend using only the MPI scheduler (mpirun) instead of utilizing multithreading (presumably via BLAS or OpenMP?), which makes sense. However, my hunch is that it's not as straightforward as simply calling mpirun on a Python script (which I tried and failed to see any speedup). Do we need to add our own MPI code, or is most of the distributed processing built into the Firedrake routines?

Error when running Lid Driven Cavity demo

Hello!
This looks like a very interesting project. I am also working on the Navier-Stokes equations with FEniCS :)

I cloned the repo and tried to run the lid driven cavity demo as it is, but it does not seem to converge. Am I doing something wrong?

(fenics_env) stefano@r11e:~/git/VarMINT$ python demos/regularizedLDC.py 
Calling FFC just-in-time (JIT) compiler, this may take some time.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 5.077e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 4.238e-02 (tol = 1.000e-10) r (rel) = 8.346e-04 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 5.532e-03 (tol = 1.000e-10) r (rel) = 1.089e-04 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 4.750e-05 (tol = 1.000e-10) r (rel) = 9.355e-07 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 5: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 6: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 7: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 8: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 9: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 10: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 11: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 12: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 13: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 14: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 15: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 16: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 17: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 18: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 19: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 20: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 21: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 22: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 23: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 24: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 25: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 26: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 27: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 28: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 29: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 30: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 31: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 32: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 33: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 34: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 35: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 36: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 37: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 38: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 39: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 40: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 41: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 42: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 43: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 44: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 45: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 46: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 47: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 48: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Traceback (most recent call last):
  File "demos/regularizedLDC.py", line 96, in <module>
    solve(F==0,up)
  File "/home/stot2794/anaconda3/envs/fenics_env/lib/python3.8/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/stot2794/anaconda3/envs/fenics_env/lib/python3.8/site-packages/dolfin/fem/solving.py", line 266, in _solve_varproblem
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     [email protected]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  c5b9b269f4a6455a739109e3a66e036b5b8412f5
*** -------------------------------------------------------------------------

Bilinear Interpolation (specifically non NN interp) from Gridded Data Sets

icepack/interpolate.py

icepack/icepack/interpolate.py

Line 53 in bda151d

q.dat.data[:] = np.fromiter(f.sample(X, indexes=1),

This uses the rasterio NN interpolation. Would be good to have at least a bilinear (or other higher-order) interpolation option.

Passing rate factor calculated from temperature back into 2D flow solver

I'm trying to implement the HeatTransport3D module to update the fluidity based on ice temperature in the 2D ice shelf example. Since the heat transport module requires a 3D mesh, I've been using the extruded mesh to calculate temperature and then trying to depth average and pass the depth averaged temperature into icepack.rate_factor to update A. However, I'm getting an error related to (I believe) the data type when passing the updated fluidity into the flow solver.

Code:

flow_model = icepack.models.IceShelf()
flow_solver = icepack.solvers.FlowSolver(flow_model, dirichlet_ids=[1])

T = firedrake.Constant(255.15)
A = icepack.rate_factor(T)

h = h0.copy(deepcopy=True)
u = flow_solver.diagnostic_solve(
    velocity=u0,
    thickness=h,
    fluidity=A,
)

mesh3d = firedrake.ExtrudedMesh(mesh, layers=1)
Q0 = firedrake.FunctionSpace(mesh3d, family='CG', degree=2, vfamily='DG', vdegree=0)
V0 = firedrake.VectorFunctionSpace(mesh3d, dim=2, family='CG', degree=2, vfamily='GL', vdegree=0)
Q2 = firedrake.FunctionSpace(mesh3d, family='CG', degree=2, vfamily='GL', vdegree=4)

heat_transport = icepack.models.HeatTransport3D()
heat_solver = icepack.solvers.HeatTransportSolver(heat_transport)

E_expr = heat_transport.energy_density(firedrake.Constant(253.0), firedrake.Constant(0.0))
E = firedrake.interpolate(E_expr, Q2)
E_init = E.copy(deepcopy=True)

final_time = 400.
num_timesteps = 200
dt = final_time / num_timesteps
a = firedrake.Constant(0.0)

W = firedrake.FunctionSpace(mesh3d, family='DG', degree=2 - 1, vfamily='GL', vdegree=1)
x, y, ζ = firedrake.SpatialCoordinate(mesh3d)

for step in tqdm.trange(num_timesteps):
    h = flow_solver.prognostic_solve(
        dt,
        thickness=h,
        velocity=u,
        accumulation=a,
        thickness_inflow=h0,
    )

    T0 = heat_transport.temperature(E)
    T = firedrake.Function(Q2)
    T.project(T0)
    T2 = icepack.depth_average(T)
    A = firedrake.project(icepack.rate_factor(T2),Q)

    u = flow_solver.diagnostic_solve(
        velocity=u,
        thickness=h,
        fluidity=A,
    )

    # elevate to 3D
    U = u.dat.data_ro[:]
    S = s.dat.data_ro[:]
    H = h.dat.data_ro[:]

    u2 = firedrake.Function(V0)
    s2 = firedrake.Function(Q0)
    h2 = firedrake.Function(Q0)

    u2.dat.data[:] = U
    s2.dat.data[:] = S
    h2.dat.data[:] = H

    # calculate vertical velocity
    ω_expr = -(u2[0].dx(0) + u2[1].dx(1)) / h2 * ζ
    ω = firedrake.project(ω_expr, W)
    w = firedrake.interpolate(h2 * ω, W)

    # calculate heating rate
    ε_x, ε_z = horizontal_strain_rate(velocity=u2, surface=s2, thickness=h2), vertical_strain_rate(velocity=u2,
                                                                                                   thickness=h2)
    τ_x, τ_z = stresses(strain_rate_x=ε_x, strain_rate_z=ε_z, fluidity=A)
    q = inner(τ_x, ε_x) + inner(τ_z, ε_z)
    qm = firedrake.interpolate(q, Q2)

    # evolve energy density
    E = heat_solver.solve(
        dt,
        energy=E,
        velocity=u2,
        vertical_velocity=w,
        thickness=h2,
        energy_inflow=E_init,
        energy_surface=E_init,
        heat=qm,
        heat_bed=firedrake.Constant(0.0),
    )

Error:

Traceback (most recent call last):
File "/firedrake/src/firedrake/firedrake/function.py", line 452, in float
raise ValueError("Can only cast scalar 'Real' Functions to float.")
ValueError: Can only cast scalar 'Real' Functions to float.
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "/firedrake/src/PyOP2/pyop2/utils.py", line 237, in verify_reshape
a = np.asarray(data, dtype=t)
ValueError: setting an array element with a sequence.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/firedrake/src/firedrake/firedrake/constant.py", line 148, in assign
self.dat.data = value
File "/firedrake/src/PyOP2/pyop2/types/glob.py", line 136, in data
self._data[:] = utils.verify_reshape(value, self.dtype, self.dim)
File "/firedrake/src/PyOP2/pyop2/utils.py", line 239, in verify_reshape
raise DataValueError("Invalid data: cannot convert to %s!" % dtype)
pyop2.exceptions.DataValueError: Invalid data: cannot convert to float64!

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "/firedrake/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3505, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "", line 1, in
runfile('/CDM_Model/230525_ImplementHeatTransport_IceShelf_2D_TwoWayCoupling.py', wdir='/CDM_Model')
File "/Applications/PyCharm.app/Contents/plugins/python/helpers/pydev/_pydev_bundle/pydev_umd.py", line 198, in runfile
pydev_imports.execfile(filename, global_vars, local_vars) # execute the script
File "/Applications/PyCharm.app/Contents/plugins/python/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "/CDM_Model/230525_ImplementHeatTransport_IceShelf_2D_TwoWayCoupling.py", line 153, in
u = flow_solver.diagnostic_solve(
File "/firedrake/src/icepack/icepack/solvers/flow_solver.py", line 138, in diagnostic_solve
return self._diagnostic_solver.solve(**kwargs)
File "/firedrake/src/icepack/icepack/solvers/flow_solver.py", line 213, in solve
self._fields[name].assign(field)
File "petsc4py/PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/firedrake/src/firedrake/firedrake/adjoint/constant.py", line 50, in wrapper
ret = assign(self, *args, **kwargs)
File "/firedrake/src/firedrake/firedrake/constant.py", line 151, in assign
raise ValueError(e)
ValueError: Invalid data: cannot convert to float64!

Managing data interpolation onto extruded meshes

Hello!

We have been developing an XZ and an XYZ hybrid model for Petermann Glacier in Northern Greenland, and as we set-up our initial conditions, we are having difficult interpolating observations onto FunctionSpaces and VectorFunctionSpaces that exist on extruded meshes.

There appears to be a function in the icepack utilities (lift3D) designed to handle a conversion from unstructured meshes to their extruded counterpart. Ideally, we could interpolate data onto the unstructured mesh, run a "lift" style function (maybe named something different to accommodate 1D to 2D lifting as well as 2D to 3D lifting) that takes the interpolated FunctionSpace on the unstructured mesh and the extruded mesh and pops out a new FunctionSpace on the extruded mesh?

Then, something that would be great is a tool to modify the layers on our new extruded mesh FunctionSpace. I think right now, the lift3d function doesn't handle multi-layer extrusion properly, so something to think about.

Thanks for any help you can provide!

Nick

change where measures (dx, ds) are applied

Right now, every model class (IceShelf, IceStream, etc.) has methods to calculate parts of the action for the diagnostic physics. For example several models call out to icepack.models.viscosity.viscosity_depth_averaged to calculate the viscous part of the action.

Each part of the action is either an integral over either the entire domain, in which case the integrand gets multiplied by the measure dx, or over the boundary, in which case the integrand gets multiplied by the measure ds. The method action creates a symbolic representation of these parts, adds them altogether, and returns the sum.

Instead, each part should return the integrand, and the action method should multiply by the right measure. The reason for this is that it becomes much easier to calculate, say, the internal viscous stress as the derivative of this part of the action. The use case for this change is to make the damage mechanics solver that @hoffmaao wrote more flexible. This change will make the process the same even when users substitute their own parameterization for the viscosity or other components, which isn't the case right now (see here). There are probably other physics models (e.g. a future calving model) that need to know the stress.

This change wouldn't actually alter what's being computed, just where the work happens. There will be a breaking change in the API for the damage solver, which would need to be initialized with the viscosity function. @hoffmaao any thoughts?

update to using new version of pygmsh

Version 7+ of pygmsh has built-in support for local mesh refinement, among other shiny new things. I pinned icepack to use version 6 or less in order to avoid breaking the API. This should be updated to support version 7, possibly with some conditionals so everything still works with version 6.

flush needed to output "ICEPACK COMPLETE" on OSX gnu compiler

When I compiled with netcdf, I had to add a
flush(nu_diag)
at the end of the program icedrv, otherwise I didn't get the last output written to the log file, including the ICEPACK COMPLETED statement, the stdio would echo "ICEPACK run did not complete".

I tried making it a
call icepack_warnings_flush(nu_diag)
but it didn't work

I'm using my macbook pro with Monterey OSX and gnu compiler. I didn't have this problem until I compiled with netcdf.
Otherwise all is well getting it to run on my laptop. Thanks, Cecilia Bitz

more specific install/ update/ reinstall instructions (esp. for editors)

In the documentation, it would be helpful to have install instructions for scenarios that are relevant to those who edit icepack. For example, someone who has been editing icepack needs to reinstall firedrake, but they don't want to reinstall icepack. What do they do?

Interpolating raster data onto 3D mesh

To ease switching between SSA and Hybrid models, it would be nice to interpolate from a 2D raster to the horizontal coordinates of a 3D mesh. Currently, using icepack.interpolate on a 3D mesh raises

ValueError: The requested sample points xi have dimension 3, but this RegularGridInterpolator has dimension 2

Since datasets like surface and bed have no depth-variability, we should just be able to interpolate using the x and y coordinates of the mesh. I guess depth-invariance is already assumed with vfamily="R", vdegree=0 on the FunctionSpace anyway?

I am happy to deal with a solution (take the first two coordinates of of the mesh), but wanted to check first--I guess this could violate some design principle in that it assumes a coordinate ordering.

add tools for mesh adaptation

Adapting the numerical grid to a particular solution is a really common pain point. There is ongoing development to add a Firedrake interface to the library Pragmatic but the timeline for this is unclear. Other packages include interfaces to tools like YAMS. We should offer something that will help users do this without lots of manual tuning.

We could implement an r-refinement scheme entirely in Firedrake. This involves solving some PDE for the mesh coordinates. For an introduction to r-adaptivity, see this book or this review. One promising approach is to formulate mesh adaptation through optimal transport, which leads to a fully nonlinear elliptic PDE called the Monge-Ampère equation. Chris Budd and coauthors have done a lot of work on optimal transport for mesh adaptation; see these slides for an introduction. Colin Cotter wrote a demo for the Firedrake documentation to solve the Monge-Ampère equation. Solving the MA equation numerically is challenging but Budd and others have found experimentally that you can do a pretty bad job in an absolute sense while still getting a relatively good numerical grid.

Whether we use optimal transport or another r-refinement scheme, we'll also have to decide on what the appropriate monitor function is. A common choice is the arc-length type functional which, in our case, would be adapting to either high strain rates, surface slopes, or both. Another obvious candidate is a monitor function based on distance to the grounding line.

Running list of useful-looking papers:

• https://doi.org/10.1007/s10092-014-0127-7
• https://doi.org/10.1016/j.jcp.2012.11.014
• https://doi.org/10.1016/j.jcp.2010.10.020
• https://doi.org/10.1137/16M1109515

rate_factor calculated using temperature output of HeatTransportSolver

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)

Code:

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)
)

error:

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!

import fails without ROL

Flagged by the Firedrake CI:

https://github.com/firedrakeproject/firedrake/runs/5275898058?check_suite_focus=true

I think ROL just needs to be added as an Icepack dependency, or have a conditional import of ROL componenets

Functions for surface/bed value extraction

Perhaps I am just missing it, but is there an easy way to extract values on the surface or bed that is compatible with Firedrake adjoint? I'd like to do inversions with HybridModel, and so need to define the loss functional just at the surface. As of now, if I just do this the naive way and define the loss functional by subtracting observations and model, I compare the depth-averaged velocity to the observations. I tried extraction using the following

def extract_surface(q_in):
    mesh_x = q_in.ufl_domain()._base_mesh
    shape = q_in.ufl_shape
    if len(shape) == 0:
        VLin = firedrake.FunctionSpace(
            q_in.ufl_domain(), "CG", 2, vfamily="CG", vdegree=1
        )
        q_targ = firedrake.interpolate(q_in, VLin)
        element_xz = q_targ.ufl_element()
        element_x = element_xz.sub_elements()[0]
    else:
        VLin = firedrake.VectorFunctionSpace(
            q_in.ufl_domain(), "CG", 2, dim=shape[0], vfamily="CG", vdegree=1
        )
        q_targ = firedrake.interpolate(q_in, VLin)
        element_xz = q_targ.ufl_element()
        element_xy = element_xz.sub_elements()[0].sub_elements()[0]
        element_x = firedrake.VectorElement(element_xy, dim=shape[0])

    Q_x = firedrake.FunctionSpace(mesh_x, element_x)
    q_x = firedrake.Function(Q_x)
    if len(shape) > 0:
        q_x.dat.data[:] = q_targ.dat.data[1::2, :]
    else:
        q_x.dat.data[:] = q_targ.dat.data[1::2]
    return q_x

but this is breaking the inversions (I think since overriding the data screws up the adjoint tape).

I also tried defining a 2d mesh in 3d space and interpolating from the 3d function space, but that produced errors around needing a VertexOnlyMesh.

Is there a solution here that I am missing, or a suggestion on how to proceed?

pip3 fails to install ROL

Hi,

I’m trying to install ROL so I can use the data assimilation capabilities of icepack but I haven't being able to install this even from source. I tried different ways and they all fail

pip3 install ROL==0.0.16

pip3 install ROL

git clone [email protected]:pyrol/pyrol.git
cd pyrol
git submodule update --init
pip3 install -e .

I also tried to build this manually via

git clone [email protected]:pyrol/pyrol.git
$ cd pyrol
$ python setup.py install

I don't have a problem with roltrilinos, I have Ubuntu 20.04.3 the latest version of pip3, CMake (3.27.4), but I still get very weird installation errors. The log is very long so I printed into a LOG_FILE & attached here (sorry about that). If anyone has another way of installing ROL, and link it to the firedrake virtualenv I would really appreciate the help. Thanks!

LOG_FILE.txt