age model test updates

.circleci/config.yml

@@ -24,15 +24,15 @@ commands:
 jobs:
   python36:
     docker:
-      - image: icepack/firedrake-python3.6:0.5.3
+      - image: icepack/firedrake-python3.6:0.5.5
     working_directory: ~/icepack
     steps:
       - build
       - unittest
       - realtest
   python38:
     docker:
-      - image: icepack/firedrake-python3.8:0.5.3
+      - image: icepack/firedrake-python3.8:0.5.5
     working_directory: ~/icepack
     steps:
       - build

README.md

@@ -1,5 +1,5 @@
 
-# icepack [![CircleCI](https://circleci.com/gh/icepack/icepack/tree/master.svg?style=svg)](https://circleci.com/gh/icepack/icepack/tree/master) [![codecov](https://codecov.io/gh/icepack/icepack/branch/master/graph/badge.svg)](https://codecov.io/gh/icepack/icepack)
+# icepack [![CircleCI](https://circleci.com/gh/icepack/icepack/tree/master.svg?style=svg)](https://circleci.com/gh/icepack/icepack/tree/master) [![codecov](https://codecov.io/gh/icepack/icepack/branch/master/graph/badge.svg)](https://codecov.io/gh/icepack/icepack) [![DOI](https://zenodo.org/badge/45697304.svg)](https://zenodo.org/badge/latestdoi/45697304)
 
 icepack is a library for modeling the flow of ice sheets and glaciers using the finite element method.
 For more information and installation instructions, see the [project webpage](https://icepack.github.io).

icepack/datasets.py

@@ -135,7 +135,7 @@ def fetch_larsen_outline():
 def fetch_mosaic_of_antarctica():
     r"""Fetch the MODIS optical image mosaic of Antarctica"""
     return nsidc_data.fetch(
-        "moa750_2009_hp1_v01.1.tif.gz",
+        "moa750_2009_hp1_v02.0.tif.gz",
         downloader=_earthdata_downloader,
         processor=pooch.Decompress(),
     )

icepack/models/__init__.py

@@ -16,6 +16,7 @@
 from icepack.models.damage_transport import DamageTransport
 from icepack.models.heat_transport import HeatTransport3D
 from icepack.models.shallow_ice import ShallowIce
+from icepack.models.age_transport import AgeTransport
 
 __all__ = [
     "IceShelf",

icepack/models/age_transport.py

@@ -0,0 +1,76 @@
+# Copyright (C) 2020-2021 by Andrew Hoffman <hoffmaao@uw.edu>
+#
+# This file is part of icepack.
+#
+# icepack is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The full text of the license can be found in the file LICENSE in the
+# icepack source directory or at <http://www.gnu.org/licenses/>.
+
+r"""Description of the age transport model.
+"""
+
+from operator import itemgetter
+import firedrake
+from firedrake import (
+    inner, grad, div, dx, ds, ds_tb, ds_v, dS_v, dS_h, dS, sqrt, det, min_value, max_value, conditional
+)
+from icepack.constants import year
+from icepack.utilities import eigenvalues, vertical_velocity, add_kwarg_wrapper
+
+def velocity_3D(**kwargs):
+    u, h = itemgetter("velocity", "thickness")(kwargs)
+
+    Q = h.function_space()
+    V = u.function_space()
+    mesh = Q.mesh()
+    xdegree_u, zdegree_u = u.ufl_element().degree()
+    W = firedrake.FunctionSpace(mesh, family='CG', degree=2, vfamily='GL', vdegree=1+V.ufl_element().degree()[1])
+    w = firedrake.interpolate(vertical_velocity(u,h),W)
+    V3D = firedrake.VectorFunctionSpace(mesh, dim=3, family='CG', degree=xdegree_u, vfamily='GL',vdegree=zdegree_u+1)
+    u3D = firedrake.Function(V3D).interpolate(firedrake.as_vector((u[0],u[1],w)))
+    return u3D 
+
+
+class AgeTransport:
+    def __init__(
+        self,
+        max_age=2000000.0,
+        velocity_3D=velocity_3D
+    ):
+        self.max_age = max_age
+        self.velocity_3D = add_kwarg_wrapper(velocity_3D)
+
+
+
+    def flux(self, **kwargs):
+        u, h, q = itemgetter("velocity", "thickness", "age")(kwargs)
+
+        u3D = self.velocity_3D(**kwargs)
+
+        Q = h.function_space()
+        U = q.function_space()
+        mesh = Q.mesh()
+
+        φ = firedrake.TestFunction(U)
+
+        xdegree_u, zdegree_u = u.ufl_element().degree()
+        degree_h = h.ufl_element().degree()[0]
+        degree = (xdegree_u + degree_h, 2 * zdegree_u + 1)
+        metadata = {'quadrature_degree': degree}
+        ice_front_ids = tuple(kwargs.pop('ice_front_ids', ()))
+
+        u3D = self.velocity_3D(**kwargs)
+
+        n = firedrake.FacetNormal(mesh)
+        u_n = max_value(0, inner(u3D, n))
+        f = q * u_n
+        q_inflow=firedrake.Constant(0.0)
+        flux_faces_v = (f('+') - f('-')) * (φ('+') - φ('-')) * dS_v
+        flux_cells = -q * div(u3D * φ) * dx
+        flux_v = q * max_value(0, inner(u3D, n)) * φ * ds_v
+        flux_tb = q * max_value(0, inner(u3D, n)) * φ * ds_tb
+        return flux_cells + flux_faces_v + flux_tb + flux_v

icepack/models/hybrid.py

@@ -217,22 +217,25 @@ def action(self, **kwargs):
         ice_front_ids = tuple(kwargs.pop("ice_front_ids", ()))
         side_wall_ids = tuple(kwargs.pop("side_wall_ids", ()))
 
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
+        ds_b = firedrake.ds_b(domain=mesh, metadata=metadata)
+        ds_v = firedrake.ds_v(domain=mesh)
+
         viscosity = self.viscosity(**kwargs) * dx
         gravity = self.gravity(**kwargs) * dx
         friction = self.friction(**kwargs) * ds_b
 
-        ds_w = ds_v(domain=mesh, subdomain_id=side_wall_ids)
-        side_friction = self.side_friction(**kwargs) * ds_w
+        side_friction = self.side_friction(**kwargs) * ds_v(side_wall_ids)
         if get_mesh_axes(mesh) == "xyz":
-            penalty = self.penalty(**kwargs) * ds_w
+            penalty = self.penalty(**kwargs) * ds_v(side_wall_ids)
         else:
             penalty = 0.0
 
         xdegree_u, zdegree_u = u.ufl_element().degree()
         degree_h = h.ufl_element().degree()[0]
         degree = (xdegree_u + degree_h, 2 * zdegree_u + 1)
-        metadata = {"quadrature_degree": degree}
-        ds_t = ds_v(domain=mesh, subdomain_id=ice_front_ids, metadata=metadata)
+        ds_t = firedrake.ds_v(ice_front_ids, metadata={"quadrature_degree": degree})
         terminus = self.terminus(**kwargs) * ds_t
 
         return viscosity + friction + side_friction - gravity - terminus + penalty
@@ -243,6 +246,11 @@ def scale(self, **kwargs):
         The positive part of the action functional is used as a dimensional
         scale to determine when to terminate an optimization algorithm.
         """
+        u = kwargs["velocity"]
+        mesh = u.ufl_domain()
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
+        ds_b = firedrake.ds_b(domain=mesh, metadata=metadata)
         return self.viscosity(**kwargs) * dx + self.friction(**kwargs) * ds_b
 
     def quadrature_degree(self, **kwargs):

icepack/models/ice_shelf.py

@@ -1,4 +1,4 @@
-# Copyright (C) 2017-2020 by Daniel Shapero <shapero@uw.edu>
+# Copyright (C) 2017-2021 by Daniel Shapero <shapero@uw.edu>
 #
 # This file is part of icepack.
 #
@@ -12,7 +12,7 @@
 
 from operator import itemgetter
 import firedrake
-from firedrake import inner, grad, dx, ds
+from firedrake import inner, grad
 from icepack.constants import ice_density as ρ_I, water_density as ρ_W, gravity as g
 from icepack.models.viscosity import viscosity_depth_averaged as viscosity
 from icepack.models.friction import side_friction, normal_flow_penalty
@@ -130,15 +130,16 @@ def action(self, **kwargs):
         ice_front_ids = tuple(kwargs.pop("ice_front_ids", ()))
         side_wall_ids = tuple(kwargs.pop("side_wall_ids", ()))
 
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
+        ds = firedrake.ds(domain=mesh, metadata=metadata)
+
         viscosity = self.viscosity(**kwargs) * dx
         gravity = self.gravity(**kwargs) * dx
 
-        ds_w = ds(domain=mesh, subdomain_id=side_wall_ids)
-        side_friction = self.side_friction(**kwargs) * ds_w
-        penalty = self.penalty(**kwargs) * ds_w
-
-        ds_t = ds(domain=mesh, subdomain_id=ice_front_ids)
-        terminus = self.terminus(**kwargs) * ds_t
+        side_friction = self.side_friction(**kwargs) * ds(side_wall_ids)
+        penalty = self.penalty(**kwargs) * ds(side_wall_ids)
+        terminus = self.terminus(**kwargs) * ds(ice_front_ids)
 
         return viscosity + side_friction - gravity - terminus + penalty
 
@@ -148,6 +149,8 @@ def scale(self, **kwargs):
         The positive part of the action functional is used as a dimensional
         scale to determine when to terminate an optimization algorithm.
         """
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
         return self.viscosity(**kwargs) * dx
 
     def quadrature_degree(self, **kwargs):

icepack/models/ice_stream.py

@@ -12,7 +12,7 @@
 
 from operator import itemgetter
 import firedrake
-from firedrake import inner, dx, ds
+from firedrake import inner
 from icepack.constants import ice_density as ρ_I, water_density as ρ_W, gravity as g
 from icepack.models.viscosity import viscosity_depth_averaged as viscosity
 from icepack.models.friction import bed_friction, side_friction, normal_flow_penalty
@@ -108,19 +108,21 @@ def action(self, **kwargs):
         ice_front_ids = tuple(kwargs.pop("ice_front_ids", ()))
         side_wall_ids = tuple(kwargs.pop("side_wall_ids", ()))
 
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
+        ds = firedrake.ds(domain=mesh, metadata=metadata)
+
         viscosity = self.viscosity(**kwargs) * dx
         friction = self.friction(**kwargs) * dx
         gravity = self.gravity(**kwargs) * dx
 
-        ds_w = ds(domain=mesh, subdomain_id=side_wall_ids)
-        side_friction = self.side_friction(**kwargs) * ds_w
+        side_friction = self.side_friction(**kwargs) * ds(side_wall_ids)
         if get_mesh_axes(mesh) == "xy":
-            penalty = self.penalty(**kwargs) * ds_w
+            penalty = self.penalty(**kwargs) * ds(side_wall_ids)
         else:
             penalty = 0.0
 
-        ds_t = ds(domain=mesh, subdomain_id=ice_front_ids)
-        terminus = self.terminus(**kwargs) * ds_t
+        terminus = self.terminus(**kwargs) * ds(ice_front_ids)
 
         return viscosity + friction + side_friction - gravity - terminus + penalty
 
@@ -130,6 +132,8 @@ def scale(self, **kwargs):
         The positive part of the action functional is used as a dimensional
         scale to determine when to terminate an optimization algorithm.
         """
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
         return (self.viscosity(**kwargs) + self.friction(**kwargs)) * dx
 
     def quadrature_degree(self, **kwargs):

icepack/models/shallow_ice.py

@@ -1,4 +1,4 @@
-# Copyright (C) 2020 by Jessica Badgeley <badgeley@uw.edu>
+# Copyright (C) 2020-2021 by Jessica Badgeley <badgeley@uw.edu>
 #
 # This file is part of icepack.
 #
@@ -12,7 +12,7 @@
 
 from operator import itemgetter
 import firedrake
-from firedrake import inner, grad, dx
+from firedrake import inner, grad
 from icepack.constants import ice_density as ρ_I, gravity as g, glen_flow_law as n
 from icepack.models.mass_transport import Continuity
 from icepack.utilities import add_kwarg_wrapper
@@ -129,6 +129,8 @@ def action(self, **kwargs):
             All other keyword arguments will be passed on to the
             'mass', 'gravity' and 'penalty' functionals
         """
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
         mass = self.mass(**kwargs) * dx
         gravity = self.gravity(**kwargs) * dx
         penalty = self.penalty(**kwargs) * dx
@@ -140,13 +142,14 @@ def scale(self, **kwargs):
         The positive part of the action functional is used as a dimensional
         scale to determine when to terminate an optimization algorithm.
         """
+        metadata = {"quadrature_degree": self.quadrature_degree(**kwargs)}
+        dx = firedrake.dx(metadata=metadata)
         return (self.mass(**kwargs) + self.penalty(**kwargs)) * dx
 
     def quadrature_degree(self, **kwargs):
         r"""Return the quadrature degree necessary to integrate the action
         functional accurately"""
-        u, h, s = itemgetter("velocity", "thickness", "surface")(kwargs)
+        u, h = itemgetter("velocity", "thickness")(kwargs)
         degree_u = u.ufl_element().degree()
         degree_h = h.ufl_element().degree()
-        degree_s = s.ufl_element().degree()
-        return int((n + 1) * degree_h + n * degree_s + degree_u)
+        return int((2 * n + 1) * degree_h + degree_u)

icepack/registry-nsidc.txt

@@ -1,6 +1,6 @@
 antarctic_ice_vel_phase_map_v01.nc fa0957618b8bd98099f4a419d7dc0e3a2c562d89e9791b4d0ed55e6017f52416 https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0754.001/1996.01.01/antarctic_ice_vel_phase_map_v01.nc
 BedMachineAntarctica_2020-07-15_v02.nc 71ff97f8fb034fbfe235cd1cad044eeea2ae1543e82115afbd886a24c1e71069 https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0756.002/1970.01.01/BedMachineAntarctica_2020-07-15_v02.nc
-moa750_2009_hp1_v01.1.tif.gz 90d1718ea0971795ec102482c47f308ba08ba2b88383facb9fe210877e80282c https://daacdata.apps.nsidc.org/pub/DATASETS/nsidc0593_moa2009/geotiff/moa750_2009_hp1_v01.1.tif.gz
+moa750_2009_hp1_v02.0.tif.gz 90d1718ea0971795ec102482c47f308ba08ba2b88383facb9fe210877e80282c https://daacdata.apps.nsidc.org/pub/DATASETS/nsidc0593_moa2009_v02/geotiff/moa750_2009_hp1_v02.0.tif.gz
 
 greenland_vel_mosaic200_2015_2016_vx_v02.1.tif 77b8eb65a4718da055bb048b75c35b48ca43d76b5bffb650932128d60ed28598 https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0478.002/2015.09.01/greenland_vel_mosaic200_2015_2016_vx_v02.1.tif
 greenland_vel_mosaic200_2015_2016_vy_v02.1.tif fb5dbc07d032de9b1bdb0b990ed02a384964d73a150529515038139efb1e3193 https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0478.002/2015.09.01/greenland_vel_mosaic200_2015_2016_vy_v02.1.tif

icepack/solvers/__init__.py

@@ -13,11 +13,13 @@
 from .flow_solver import FlowSolver, ImplicitEuler, LaxWendroff
 from .heat_transport import HeatTransportSolver
 from .damage_solver import DamageSolver
+from .age_solver import AgeSolver
 
 __all__ = [
-    "FlowSolver",
-    "ImplicitEuler",
-    "LaxWendroff",
-    "HeatTransportSolver",
-    "DamageSolver",
+    'FlowSolver',
+    'ImplicitEuler',
+    'LaxWendroff',
+    'HeatTransportSolver',
+    'DamageSolver',
+    'AgeSolver'
 ]