released 2026.3.13 with fixed firedrake version compatibility issues

Author: nicholasmr

src/Makefile

@@ -80,10 +80,10 @@ devinstall: specfabpy
 
 # Firedrake/icepack docker enviroment
 icepack: specfabpy
-	@echo "*** You can also just run: pip3 install specfabpy ***"
-	#pip3 install setuptools==61 cartopy cmasher
-	#pip3 install . 
+	@echo "*** Alternatively: pip3 install specfabpy ***"
+	pip3 install cmasher tabulate termcolor cartopy scipy
 	python3 setup.py install
+	#pip3 install . 
 
 specfabpy: $(SPECFAB).o
 	@echo "*** Archiving objects into $(PWD)/libspecfab.a";

src/specfabpy/fenics/steadyCPO.py

@@ -7,6 +7,7 @@
 
 import os, sys, copy, code # code.interact(local=locals())
 import numpy as np
+from scipy.spatial.transform import Rotation as R
 import xarray, pickle, pyproj # pyproj needed!
 from tabulate import tabulate
 
@@ -195,7 +196,7 @@ def get_mesh(self, mapscale=1):
         cells = mesh.cells()
         return (mesh, boundaries, Q,Q2,V, coords,cells)
 
-    def triang(self, coords, cells, mapscale=1):
+    def get_triang(self, coords, cells, mapscale=1):
 
         return tri.Triangulation(mapscale*coords[0], mapscale*coords[1], triangles=cells)
 
@@ -316,7 +317,28 @@ def set_solution(self, probname):
     def set_inputs(self):
 
         self.coords,self.cells, self.ux,self.uy,self.umag,self.epsE, self.S,self.B,self.H,self.mask = self.npinputs()
-        self.triang = self.triang(self.coords, self.cells, mapscale=self.mapscale)
+        self.triang = self.get_triang(self.coords, self.cells, mapscale=self.mapscale)
+
+    def set_adjustedcoords(self, xoffset=0, yoffset=0, rotang=0):
+    
+        ### Translate points so that (x0,y0) is new origin
+        
+        xt, yt = self.coords[0]-xoffset, self.coords[1]-yoffset # translated coords
+        x0t, y0t = self.x0-xoffset, self.y0-yoffset # mesh corners
+        x1t, y1t = self.x1-xoffset, self.y1-yoffset # mesh corners
+
+        # rotate points around new origin
+        if rotang:
+            Q = R.from_euler('z', rotang, degrees=True).as_matrix()[0:2,0:2]
+            xf, yf = np.array([ np.matmul(Q,[xt[ii],yt[ii]]) for ii in range(len(xt[:])) ]).T
+            self.ux, self.uy = np.array([ np.matmul(Q,[self.ux[ii],self.uy[ii]]) for ii in range(len(self.ux[:])) ]).T
+        else:
+            xf, yf = xt, yt
+    
+        self.coords[0,:], self.coords[1,:] = xf, yf
+        self.triang = self.get_triang(self.coords, self.cells, mapscale=self.mapscale)
+
+        return Q
 
     def bmesh(self, mapscale=1):
 
@@ -386,6 +408,7 @@ def plot_generic(self, ax, F, kw_tcf=dict(), kw_cb=dict(), kw_cax=dict()):
             kw_tcf.pop('vcenter')
         cs = ax.tricontourf(self.triang, F, **kw_tcf)
         hcb = plt.colorbar(cs, cax=self.newcax(ax, **kw_cax), **kw_cb)
+        hcb.ax.yaxis.label.set_multialignment('center')
         return (cs, hcb)
 
     def plot_velocities(self, ax, kw_cb=dict(label=r'$u$ (m/yr)'), kw_cax=dict(), \
@@ -461,8 +484,8 @@ def savefig(self, fig, fname, dpi=150, pad_inches=0.1, bbox_inches='tight'):
 
         fig.savefig(fname, dpi=dpi, pad_inches=pad_inches, bbox_inches=bbox_inches)
 
-    def setupaxis(self, ax, boundaries=False, floating=True, mesh=False, bgcolor='0.85', showyaxis=True, \
-             xlims=None, ylims=None, xticks_major=None, xticks_minor=None, yticks_major=None, yticks_minor=None):
+    def setupaxis(self, ax, boundaries=False, floating=True, mesh=False, bgcolor='0.85', showxaxis=True, showyaxis=True, axform='square', \
+             xlims=None, ylims=None, xticks_major=None, xticks_minor=None, yticks_major=None, yticks_minor=None, ylabelpad=None):
 
         legh, legt = [], []
 
@@ -488,16 +511,19 @@ def setupaxis(self, ax, boundaries=False, floating=True, mesh=False, bgcolor='0.
             legh.append(Line2D([0], [0], color=self.c_floating, lw=2))
             legt.append('Floating')
 
-        ax.axis('square')
-        ax.set_xlabel(r'$x$ (km)')
+        if axform is not None: ax.axis(axform)
 
-        if xticks_major is not None: ax.set_xticks(xticks_major)
-        if xticks_minor is not None: ax.set_xticks(xticks_minor, minor=True)
-            
+        if showxaxis:
+            if xticks_major is not None: ax.set_xticks(xticks_major)
+            if xticks_minor is not None: ax.set_xticks(xticks_minor, minor=True)
+            ax.set_xlabel(r'$x$ (km)')
+        else:
+            ax.tick_params('x', labelbottom=False)
+                        
         if showyaxis:
             if yticks_major is not None: ax.set_yticks(yticks_major)
             if yticks_minor is not None: ax.set_yticks(yticks_minor, minor=True)
-            ax.set_ylabel(r'$y$ (km)')
+            ax.set_ylabel(r'$y$ (km)', labelpad=ylabelpad)
         else:
             ax.tick_params('y', labelleft=False)
 

src/specfabpy/firedrake/ice.py

@@ -189,9 +189,9 @@ def weakform(self, s, u,S, dt, iota=None, Gamma0=None, Lambda0=None, zeta=0):
         ENABLE_REG  = True
 
         # Flattened strain-rate and spin tensors for accessing them per node
-        Df = project( sym(grad(u)), self.G).vector()[:] # strain rate
-        Wf = project(skew(grad(u)), self.G).vector()[:] # spin
-        Sf = project(S, self.G).vector()[:] # deviatoric stress
+        Df = project( sym(grad(u)), self.G).dat.data_ro # strain rate
+        Wf = project(skew(grad(u)), self.G).dat.data_ro # spin
+        Sf = project(S, self.G).dat.data_ro # deviatoric stress
 
         # Same but in 3D for fabric problem
         D3 = self.mat3d(Df) # [node,3,3]
@@ -205,16 +205,16 @@ def weakform(self, s, u,S, dt, iota=None, Gamma0=None, Lambda0=None, zeta=0):
 
         # Populate entries of dynamical matrices *row-wise* (row index is ii)
         for ii in self.srrng:
-            if ENABLE_LROT: self.Mrr_LROT[ii].vector()[:]     = M_LROT[:,ii,:] # all nodes, row=ii, all columns
-            if ENABLE_DDRX: self.Mrr_DDRX_src[ii].vector()[:] = M_DDRX_src[:,ii,:]
-            if ENABLE_REG:  self.Mrr_REG[ii].vector()[:]      = M_REG[:,ii,:]
+            if ENABLE_LROT: self.Mrr_LROT[ii].dat.data[:]     = M_LROT[:,ii,:] # all nodes, row=ii, all columns
+            if ENABLE_DDRX: self.Mrr_DDRX_src[ii].dat.data[:] = M_DDRX_src[:,ii,:]
+            if ENABLE_REG:  self.Mrr_REG[ii].dat.data[:]      = M_REG[:,ii,:]
 
         ### Construct weak form
 
         # dummy zero term to make rhs(F) work when solving steady-state problem 
         # this can probably be removed once the SSA source/sink terms are added
         s_null = Function(self.S)
-        s_null.vector()[:] = 0.0
+        s_null.dat.data[:] = 0.0
         F = dot(s_null, self.w)*dx
 
         # Time derivative
@@ -277,7 +277,7 @@ def _setaux(self, u=None):
         Set auxiliary fields
         """
         sp = project(self.s, self.Sd)
-        self.rnlm = np.array([sp.sub(ii).vector()[:] + 0j for ii in self.srrng]) # reduced form (rnlm) per node
+        self.rnlm = np.array([sp.sub(ii).dat.data_ro + 0j for ii in self.srrng]) # reduced form (rnlm) per node
         self.nlm  = np.array([self.sf.rnlm_to_nlm(self.rnlm[:,nn], self.nlm_len) for nn in self.dofs0]) # full form (nlm) per node (nlm[node,coef])
 
         if self.setextra:
@@ -302,26 +302,26 @@ def get_pfJ(self, *args, **kwargs):
         Pole figure J (pfJ) index
         """
         pfJ = Function(self.Rd)
-        pfJ.vector()[:] = sfcom.pfJ(self.nlm, *args, **kwargs)[:]
+        pfJ.dat.data[:] = sfcom.pfJ(self.nlm, *args, **kwargs)[:]
         return pfJ
 
     def get_E_CAFFE(self, u):
         """
         CAFFE model (Placidi et al., 2010)
         """
-        Df = project(sym(grad(u)), self.Gd).vector()[:]
+        Df = project(sym(grad(u)), self.Gd).dat.data_ro
         E_CAFFE = Function(self.Rd)
-        E_CAFFE.vector()[:] = self.sf.E_CAFFE_arr(self.nlm, self.mat3d(Df), *self.CAFFE_params)
+        E_CAFFE.dat.data[:] = self.sf.E_CAFFE_arr(self.nlm, self.mat3d(Df), *self.CAFFE_params)
         return E_CAFFE
 
     def get_E_EIE(self, u):
         """
         EIE model (Rathmann et al., in prep)
         *** Not yet implemented ***
         """
-#        Df = project(sym(grad(u)), self.Gd).vector()[:]
+#        Df = project(sym(grad(u)), self.Gd).dat.data_ro
         E_EIE = Function(self.Rd)
-#        E_EIE.vector()[:] = self.sf.E_EIE_arr(...)
+#        E_EIE.dat.data[:] = self.sf.E_EIE_arr(...)
         return E_EIE
 
     def get_Eij(self, ei=(), xyz_sort=False):   
@@ -348,12 +348,12 @@ def get_Eij(self, ei=(), xyz_sort=False):
         Eij_fs  = [Function(self.Rd) for _ in range(6)] # Eij enhancement tensor
 
         for ii in range(3): # mi=m1,m2,m3
-            lami_fs[ii].vector()[:] = lami[:,ii]
+            lami_fs[ii].dat.data[:] = lami[:,ii]
             for jj in range(3): # j=x,y,z
-                ei_fs[ii].sub(jj).vector()[:] = ei[ii][:,jj]
+                ei_fs[ii].sub(jj).dat.data[:] = ei[ii][:,jj]
 
         for kk in range(6): 
-            Eij_fs[kk].vector()[:] = Eij[:,kk]
+            Eij_fs[kk].dat.data[:] = Eij[:,kk]
 
         return (ei_fs, Eij_fs, lami_fs)
 
@@ -371,12 +371,12 @@ def get_lamxi(self, s=None):
 
         rows = np.arange(len(lami))
         lami_fs = [Function(self.Rd) for _ in range(3)] # a2 eigenvalues (lami)
-        lami_fs[0].vector()[:] = lami[rows, Ix]
-        lami_fs[1].vector()[:] = lami[rows, Iy]
-        lami_fs[2].vector()[:] = lami[rows, Iz]
+        lami_fs[0].dat.data[:] = lami[rows, Ix]
+        lami_fs[1].dat.data[:] = lami[rows, Iy]
+        lami_fs[2].dat.data[:] = lami[rows, Iz]
 
         dlamxy = Function(self.Rd)
-        dlamxy.vector()[:] = lami[rows, Ix] - lami[rows, Iy]
+        dlamxy.dat.data[:] = lami[rows, Ix] - lami[rows, Iy]
 
         return lami_fs, dlamxy
 

src/specfabpy/firedrake/mesh.py

@@ -66,7 +66,7 @@ def print_setup(self):
 
     def dt_CFL(self, u):
 
-        ux_max = abs(u.sub(0).vector()[:]).max()
+        ux_max = abs(u.sub(0).dat.data_ro).max()
         dt_CFL = 0.5*self.hx_min/ux_max
         return dt_CFL # CFL timestep size
 
@@ -195,7 +195,7 @@ def print_setup(self):
 
     def dt_CFL(self, u):
 
-        ux_max = abs(u.sub(0).vector()[:]).max()
+        ux_max = abs(u.sub(0).dat.data_ro).max()
         dt_CFL = 0.5*self.hx_min/ux_max
         return dt_CFL # CFL timestep size
 

src/specfabpy/plotting.py

@@ -164,6 +164,7 @@ def plotcoordaxes(ax, geo, axislabels='xi', color=c_dred, fontsize=None, negaxes
         elif axislabels == 'vuxi':    lbls = [r'$\vu{x}$',r'$\vu{y}$',r'$\vu{z}$']
         elif axislabels == 'vuei':    lbls = [r'$\vu{e}_1$',r'$\vu{e}_2$',r'$\vu{e}_3$']
         elif axislabels == 'vmi':     lbls = [r'$\vb{m}_1$',r'$\vb{m}_2$',r'$\vb{m}_3$']
+        elif axislabels == 'vunt':    lbls = [r'$\vu{t}$',r'$\vu{n}$',r'$\vu{z}$']
         else: raise ValueError('sfplt.plotcoordinateaxes(): Note sure what to do with passed axislabels.')
     elif isinstance(axislabels, list) and len(axislabels) == 3:
         lbls = axislabels

tests/firedrake/simpleshear.py

@@ -41,7 +41,7 @@
 iota = +1  # deck-of-cards behaviour for lattice rotation
 Tice = -15 # Ice temperature (deg. C) for DDRX rate factor
 
-ENABLE_DDRX = True
+ENABLE_DDRX = False
 
 ### Viscous anisotropy homogenization parameters
 
@@ -74,7 +74,7 @@
 tau = fd.project(fd.sym(fd.grad(u)), T) # assume driving stress is coaxial to strain-rate (in two-way coupling this should be modelled tau) 
 
 h_min = 1/nx
-v_max = abs(u.vector()[:]).max()
+v_max = abs(u.dat.data_ro).max()
 dt_CFL = 0.5*h_min/v_max
 dt = 2*dt_CFL # more aggresive time-stepping than CFL
 
@@ -83,7 +83,7 @@
 fabric = IceFabric(mesh, boundaries, L, **kw_num, **kw_VA) # initializes as isotropic fabric field
 fabric.set_isotropic_BCs((1,)) # isotropic ice incoming from left-hand boundary, remaining boundaries are free (no fabric fluxes)
 
-Gamma0 = None
+Gamma0 = Gamma0_list = None
 if ENABLE_DDRX:
     Gamma0      = fabric.Gamma0_Lilien23_EDC(u,Tice+273.15)
     Gamma0_list = [fabric.Gamma0_Lilien23_EDC(u,_+273.15) for _ in np.linspace(-40,Tice,3)] # list of DDRX rate factors used to gradually approach solution