diff --git a/src/specfabpy/fenics/CPO.py b/src/specfabpy/fenics/CPO.py index d4b25c9..0747cc4 100644 --- a/src/specfabpy/fenics/CPO.py +++ b/src/specfabpy/fenics/CPO.py @@ -21,7 +21,7 @@ def __init__(self, mesh, boundaries, L, nu_multiplier=1, nu_realspace=1e-3, mode ### Check args - if modelplane != 'xz': + if modelplane not in ['xz', 'xy']: raise ValueError('modelplane "%s" not supported, must be "xz"'%(modelplane)) ### Setup @@ -80,40 +80,41 @@ def __init__(self, mesh, boundaries, L, nu_multiplier=1, nu_realspace=1e-3, mode self.wr = TestFunction(self.S) # weight functions (real part of nlm coefs) self.sr_sub = split(self.sr) # for easy access of each subelement self.wr_sub = split(self.wr) # for easy access of each subelement - + ### Solution containers - + self.s = Function(self.S) # Current solution self.s0 = Function(self.S) # Previous solution - + ### Dynamical matrices - + self.nlm_dummy = np.zeros((self.nlm_len_full)) self.Mk_LROT = [ [Function(self.V) for ii in self.srng] for ii in range(4) ] # rr, ri, ir, ii self.Mk_DDRX_src = [ [Function(self.V) for ii in self.srng] for ii in range(4) ] self.Mk_CDRX = [ [Function(self.V) for ii in self.srng] for ii in range(4) ] self.Mk_REG = [ [Function(self.V) for ii in self.srng] for ii in range(4) ] - + ### Aux - + # Idealized states self.nlm_iso = [1/np.sqrt(4*np.pi)] + [0]*(self.nlm_len-1) # Isotropic and normalized state self.nlm_zero = [0]*(self.nlm_len) - + self.M_zero = np.zeros((self.numdofs, self.nlm_len, self.nlm_len)) + def initialize(self, sr=None, si=None): """ Initialize uniform CPO field """ - - if sr is None: sr, si = self.nlm_iso, self.nlm_zero - + if sr is None: + sr, si = self.nlm_iso, self.nlm_zero + s = Function(self.S) if self.modelplane=='xy': - assign(self.s.sub(0), project(Constant(sr), self.V)) # real part - assign(self.s.sub(1), project(Constant(si), self.V)) # imag part + assign(s.sub(0), project(Constant(sr), self.V)) # real part + assign(s.sub(1), project(Constant(si), self.V)) # imag part elif self.modelplane=='xz': - assign(self.s, project(Constant(sr), self.V)) # real part - + assign(s, project(Constant(sr), self.V)) # real part + self.set_state(s) def set_state(self, s, interp=False): if interp: @@ -121,16 +122,13 @@ def set_state(self, s, interp=False): else: self.s.assign(s) self.s0.assign(s) - def set_BCs(self, sr, si, domids, domain=None): """ Set boundary conditions - """ - + """ self.bcs = [] if domain is None: domain = self.boundaries - for ii, did in enumerate(domids): if self.modelplane=='xy': self.bcs += [DirichletBC(self.S.sub(0), sr[ii], domain, did)] # real part @@ -138,30 +136,22 @@ def set_BCs(self, sr, si, domids, domain=None): elif self.modelplane=='xz': self.bcs += [DirichletBC(self.S, sr[ii], domain, did)] # real part - def set_isotropic_BCs(self, domids, domain=None): """ Easy setting of isotropic boundary conditions """ - sr = [Constant(self.nlm_iso)] * len(domids) si = [Constant(self.nlm_zero)] * len(domids) self.set_BCs(sr, si, domids, domain=domain) - def evolve(self, u, S, dt, iota=+1, Gamma0=None, Lambda0=None, steadystate=False, disable_advection=False): """ The fabric solver, called to step fabric field forward in time by amount dt - """ - - if self.s is None: - raise ValueError('CPO state "w" not set. Did you forget to initialize the CPO field?') - + """ self.s0.assign(self.s) # current state self.s must be set F = self._weakform(u, S, dt, iota, Gamma0, Lambda0, steadystate=steadystate, disable_advection=disable_advection) solve(lhs(F)==rhs(F), self.s, self.bcs, solver_parameters={'linear_solver':'gmres', }) # fastest tested are: gmres, bicgstab, tfqmr --- note this is a non-symmetric system! - def solvesteady(self, u, S, iota=+1, Gamma0=None, Lambda0=None, disable_advection=False, \ dt=None, tol=1e-7, maxiter=400): @@ -184,13 +174,13 @@ def solvesteady(self, u, S, iota=+1, Gamma0=None, Lambda0=None, disable_advectio else: self.evolve(u, S, 1, steadystate=True, **kwargs_dyn) - def _weakform(self, u, S, dt, iota, Gamma0, Lambda0, zeta=0, steadystate=False, disable_advection=False): ENABLE_LROT = iota is not None ENABLE_DDRX = (Gamma0 is not None) and (S is not None) ENABLE_CDRX = Lambda0 is not None - + 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 @@ -201,25 +191,19 @@ def _weakform(self, u, S, dt, iota, Gamma0, Lambda0, zeta=0, steadystate=False, W3 = self.mat3d(Wf) S3 = self.mat3d(Sf) - # Dynamical matrices at each DOF - if ENABLE_LROT: M_LROT_nodal = np.array([self.sf.reduce_M(self.sf.M_LROT(self.nlm_dummy, D3[nn], W3[nn], iota, zeta), self.nlm_len) for nn in self.dofs] ) - if ENABLE_DDRX: M_DDRX_src_nodal = np.array([self.sf.reduce_M(self.sf.M_DDRX_src(self.nlm_dummy, S3[nn]), self.nlm_len) for nn in self.dofs] ) - M_REG_nodal = np.array([self.sf.reduce_M(self.sf.M_REG(self.nlm_dummy, D3[nn]), self.nlm_len) for nn in self.dofs] ) - -# @TODO IMPLEMENT -# # Dynamical matrices M_* at each node as np arrays (indexing is M[node,row,column]) -# M_LROT = np.array([self._M_reduced(self.sf.M_LROT, self.nlm_dummy, D3[nn], W3[nn], iota, zeta) for nn in self.dofs] ) if ENABLE_LROT else self.M_zero -# M_DDRX_src = np.array([self._M_reduced(self.sf.M_DDRX_src, self.nlm_dummy, S3[nn]) for nn in self.dofs] ) if ENABLE_DDRX else self.M_zero -# M_REG = np.array([self._M_reduced(self.sf.M_REG, self.nlm_dummy, D3[nn]) for nn in self.dofs] ) if ENABLE_REG else self.M_zero + # Dynamical matrices M_* at each node as np arrays (indexing is M[node,kk,row,column]) + M_LROT = np.array([self._M_reduced(self.sf.M_LROT, self.nlm_dummy, D3[nn], W3[nn], iota, zeta) for nn in self.dofs] ) if ENABLE_LROT else self.M_zero + M_DDRX_src = np.array([self._M_reduced(self.sf.M_DDRX_src, self.nlm_dummy, S3[nn]) for nn in self.dofs] ) if ENABLE_DDRX else self.M_zero + M_REG = np.array([self._M_reduced(self.sf.M_REG, self.nlm_dummy, D3[nn]) for nn in self.dofs] ) if ENABLE_REG else self.M_zero # Populate entries of dynamical matrices if self.modelplane=='xy': krng = range(4) # rr, ri, ir, ii elif self.modelplane=='xz': krng = range(1) # rr for ii in self.srng: for kk in krng: - if ENABLE_LROT: self.Mk_LROT[kk][ii].vector()[:] = M_LROT_nodal[:,kk,ii,:].flatten() - if ENABLE_DDRX: self.Mk_DDRX_src[kk][ii].vector()[:] = M_DDRX_src_nodal[:,kk,ii,:].flatten() - self.Mk_REG[kk][ii].vector()[:] = M_REG_nodal[:,kk,ii,:].flatten() + if ENABLE_LROT: self.Mk_LROT[kk][ii].vector()[:] = M_LROT[:,kk,ii,:].flatten() + if ENABLE_DDRX: self.Mk_DDRX_src[kk][ii].vector()[:] = M_DDRX_src[:,kk,ii,:].flatten() + if ENABLE_REG: self.Mk_REG[kk][ii].vector()[:] = M_REG[:,kk,ii,:].flatten() ### Construct weak form @@ -252,15 +236,19 @@ def _weakform(self, u, S, dt, iota, Gamma0, Lambda0, zeta=0, steadystate=False, raise ValueError('CPO(): DDRX evolution is not yet supported for modelplane="xy"') # Orientation space stabilization (hyper diffusion) - Mrr_REG, Mri_REG, Mir_REG, Mii_REG = self.Mk_REG # unpack for readability - F += -self.nu_multiplier * sum([ (dot(Mrr_REG[ii], self.sr) + dot(Mri_REG[ii], self.si))*self.wr_sub[ii]*dx for ii in self.srng]) # real part - F += -self.nu_multiplier * sum([ (dot(Mir_REG[ii], self.sr) + dot(Mii_REG[ii], self.si))*self.wi_sub[ii]*dx for ii in self.srng]) # imag part + if ENABLE_REG: + Mrr_REG, Mri_REG, Mir_REG, Mii_REG = self.Mk_REG # unpack for readability + F += -self.nu_multiplier * sum([ (dot(Mrr_REG[ii], self.sr) + dot(Mri_REG[ii], self.si))*self.wr_sub[ii]*dx for ii in self.srng]) # real part + F += -self.nu_multiplier * sum([ (dot(Mir_REG[ii], self.sr) + dot(Mii_REG[ii], self.si))*self.wi_sub[ii]*dx for ii in self.srng]) # imag part elif self.modelplane=='xz': - # Real space stabilization (Laplacian diffusion) - F = self.nu_realspace * inner(grad(self.sr), grad(self.wr))*dx # real part - + # 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 + F = dot(s_null, self.wr)*dx + # Time derivative if not steadystate: F += dtinv * dot( (self.sr-self.s0), self.wr)*dx # real part @@ -268,6 +256,9 @@ def _weakform(self, u, S, dt, iota, Gamma0, Lambda0, zeta=0, steadystate=False, # Real space advection if not disable_advection: F += dot(dot(u, nabla_grad(self.sr)), self.wr)*dx # real part + + # Real space stabilization (Laplacian diffusion) + F += self.nu_realspace * inner(grad(self.sr), grad(self.wr))*dx # real part # Lattice rotation if ENABLE_LROT: @@ -282,45 +273,39 @@ def _weakform(self, u, S, dt, iota, Gamma0, Lambda0, zeta=0, steadystate=False, F += F_src - F_sink # Orientation space stabilization (hyper diffusion) - Mrr_REG, *_ = self.Mk_REG # unpack for readability - F += -self.nu_multiplier * sum([ dot(Mrr_REG[ii], self.sr)*self.wr_sub[ii]*dx for ii in self.srng]) # real part + if ENABLE_REG: + Mrr_REG, *_ = self.Mk_REG # unpack for readability + F += -self.nu_multiplier * sum([ dot(Mrr_REG[ii], self.sr)*self.wr_sub[ii]*dx for ii in self.srng]) # real part return F - def _M_reduced(self, Mfun, *args, **kwargs): - return self.sf.reduce_M(Mfun(*args, **kwargs), self.nlm_len)[0] # full form -> reduced form of M_* + return self.sf.reduce_M(Mfun(*args, **kwargs), self.nlm_len) # full form -> reduced form of M_* - def mat3d(self, D2, dn=4): D2 = np.array([ D2[nn*dn:(nn+1)*dn] for nn in self.dofs ]) # to [node, D2 flat index 0 to 3] return sfcom.mat3d_arr(D2, self.modelplane, reshape=True) # [node,3,3] - - def get_nlm(self, x, y): - - """ - Extract CPO state gridpoint-wise (full form of nlm) - """ - - if self.modelplane=='xy': nlm = self.s.sub(0)(x,y) + 1j * self.s.sub(1)(x,y) - elif self.modelplane=='xz': nlm = self.s(x,y) + 0j - return self.sf.rnlm_to_nlm(nlm, self.nlm_len_full) if self.USE_REDUCED else nlm - - def eigenframe(self, *p, **kwargs): """ Extract a2 eigenvectors (mi) and eigenvalues (lami) at specified list of points p=([x1,x2,...], [y1,y2,...]) Note that you can back-construct a2 = lam1*np.outer(m1,m1) + lam2*np.outer(m2,m2) + lam3*np.outer(m3,m3) """ - xf, yf = np.array(p[0], ndmin=1), np.array(p[1], ndmin=1) # (point num, x/y value) nlm = np.array([self.get_nlm(xf[pp],yf[pp], **kwargs) for pp in range(len(xf))]) # (point num, nlm coefs) mi, lami = sfcom.eigenframe(nlm, symframe=self.symframe, modelplane=self.modelplane) return (mi, lami) - + + def get_nlm(self, *p, xz2xy=False): + """ + Extract CPO state gridpoint-wise (full form of nlm) + """ + if self.modelplane=='xy': nlm_ = self.s.sub(0)(*p) + 1j * self.s.sub(1)(*p) + elif self.modelplane=='xz': nlm_ = self.s(*p) + 0j + nlm = self.sf.rnlm_to_nlm(nlm_, self.nlm_len_full) if self.USE_REDUCED else nlm_ + return self.sf.rotate_nlm_xz2xy(nlm) if xz2xy else nlm + def apply_bounds(self): - """ Renormalize power spectrum if it exceeds that of the delta function (numerical overshoot) """ diff --git a/src/specfabpy/fenics/ice.py b/src/specfabpy/fenics/ice.py index 16f9663..a81bd46 100644 --- a/src/specfabpy/fenics/ice.py +++ b/src/specfabpy/fenics/ice.py @@ -11,53 +11,83 @@ class IceFabric(CPO): - def __init__(self, mesh, boundaries, *args, L=8, nu_realspace=1e-3, nu_multiplier=1, modelplane='xz', ds=None, nvec=None, \ - Eij_grain=(1,1), alpha=0, n_grain=1, symframe=-1, \ - Cij=sfconst.ice['elastic']['Bennett1968'], rho=sfconst.ice['density'], **kwargs): + def __init__( + self, mesh, boundaries, L, *args, + nu_multiplier=1, nu_realspace=1e-3, modelplane='xz', symframe=-1, + Eij_grain=(1,1), alpha=0, n_grain=1, CAFFE_params=(0.1, 10), n_EIE=3, + ds=None, nvec=None, setextra=True, + Cij=sfconst.ice['elastic']['Bennett1968'], rho=sfconst.ice['density'], **kwargs + ): + ### Check args + # ...done in CPO() and EnhancementFactor() + + ### Setup + + # CPO class + self.setextra = setextra # set Eij, E_CAFFE, pfJ, etc. on every state update? super().__init__(mesh, boundaries, L, nu_multiplier=nu_multiplier, nu_realspace=nu_realspace, modelplane=modelplane, ds=ds, nvec=nvec) - self.initialize(sr=None) # isotropic - self.set_BCs([], [], []) # no BCs + # specfab self.grain_params = (Eij_grain, alpha, n_grain) + self.CAFFE_params = CAFFE_params # (Emin, Emax) + self.n_EIE = n_EIE self.Lame_grain = self.sf.Cij_to_Lame_tranisotropic(Cij) self.rho = rho + # EnhancementFactor class self.symframe = symframe self.enhancementfactor = EnhancementFactor(mesh, L, symframe=symframe, modelplane=modelplane) - self.update_Eij() - - def get_state(self, *args, **kwargs): - return self.get_nlm(*args, **kwargs) # alias + + ### Finish up + + self.initialize() # isotropic + self.set_BCs([], [], []) # no BCs + + + def initialize(self, *args, **kwargs): + super().initialize(*args, **kwargs) + self._setaux() def set_state(self, *args, **kwargs): super().set_state(*args, **kwargs) - self.update_Eij() + self._setaux() - def evolve(self, *args, **kwargs): - super().evolve(*args, **kwargs) - self.update_Eij() + def evolve(self, u, *args, **kwargs): + super().evolve(u, *args, **kwargs) + self._setaux(u=u) def solvesteady(self, u, S, **kwargs): super().solvesteady(u, S, **kwargs) - self.update_Eij() + self._setaux(u=u) - def update_Eij(self): - self.mi, self.Eij, self.lami = self.enhancementfactor.Eij_tranisotropic(self.s, *self.grain_params, ei=()) - self.xi, self.Exij, _ = self.enhancementfactor.Eij_tranisotropic(self.s, *self.grain_params, ei=np.eye(3)) - self.pfJ = self.enhancementfactor.pfJ(self.s) - # ... unpack - self.m1, self.m2, self.m3 = self.mi # rheological symmetry directions - self.E11, self.E22, self.E33, self.E23, self.E31, self.E12 = self.Eij # eigenenhancements - self.Exx, self.Eyy, self.Ezz, self.Eyz, self.Exz, self.Exy = self.Exij # Cartesian enhancements - self.lam1, self.lam2, self.lam3 = self.lami # fabric eigenvalues \lambda_i (not a2 eigenvalues unless symframe=-1) - - def E_CAFFE(self, u, *args, **kwargs): - return self.enhancementfactor.E_CAFFE(self.s, u, *args, **kwargs) + def _setaux(self, u=None): + if self.setextra: + self.mi, self.Eij, self.lami = self.get_Eij(ei=()) + self.xi, self.Exij, _ = self.get_Eij(ei=np.eye(3)) + # unpack above eigenframe fields for convenience + self.m1, self.m2, self.m3 = self.mi # rheological symmetry directions + self.E11, self.E22, self.E33, self.E23, self.E31, self.E12 = self.Eij # eigenenhancements + self.Exx, self.Eyy, self.Ezz, self.Eyz, self.Exz, self.Exy = self.Exij # Cartesian enhancements + self.lam1, self.lam2, self.lam3 = self.lami # fabric eigenvalues \lambda_i (not a2 eigenvalues unless symframe=-1) + # additional fields... + self.pfJ = self.get_pfJ() + if u is not None: + self.E_CAFFE = self.get_E_CAFFE(u) + self.E_EIE = self.get_E_EIE(u, self.Eij, self.mi) - def E_EIE(self, u, *args, **kwargs): - return self.enhancementfactor.E_EIE(u, self.Eij, self.mi, *args, **kwargs) + def get_pfJ(self, *args, **kwargs): + return self.enhancementfactor.pfJ(self.s, *args, **kwargs) + def get_E_CAFFE(self, u, *args, **kwargs): + return self.enhancementfactor.E_CAFFE(self.s, u, Emin=self.CAFFE_params[0], Emax=self.CAFFE_params[1]) + + def get_E_EIE(self, u, Eij, mi, *args, **kwargs): + return self.enhancementfactor.E_EIE(u, Eij, mi, self.n_EIE) + + def get_Eij(self, ei=()): + return self.enhancementfactor.Eij_tranisotropic(self.s, *self.grain_params, ei=()) + def get_elastic_velocities(self, x,y, theta,phi, alpha=1): nlm = self.get_state(x,y) vS1, vS2, vP = sf__.Vi_elastic_tranisotropic(nlm, alpha,self.Lame_grain,self.rho, theta,phi) # calculate elastic phase velocities using specfab diff --git a/src/specfabpy/firedrake/ice.py b/src/specfabpy/firedrake/ice.py index d83da72..44d37f8 100644 --- a/src/specfabpy/firedrake/ice.py +++ b/src/specfabpy/firedrake/ice.py @@ -1,12 +1,12 @@ #!/usr/bin/python3 # Nicholas Rathmann and Daniel Shapero, 2024 -r""" +""" Firedrake all-in-one interface for ice fabric dynamics, viscous anisotropy, etc. ------- +----------- NOTICE ------- +----------- The 2D ice-flow model plane is assumed to be xz. This has the benifit of reducing the DOFs of the fabric evolution problem to involve only real numbers (real-valued state vector, s), considerably speeding up the solver. If your problem is in fact xy, nothing changes in the way this class is used; the (M)ODFs etc will simply reflect to assumption that flow is in xz rather than xy. @@ -25,44 +25,65 @@ iota : plastic spin free parameter; iota=1 => deck-of-cards behaviour Gamma0 : DDRX rate factor so that Gamma=Gamma0*(D-), where D is the deformability Lambda0 : CDRX rate factor -""" -""" -TODO: +----------- +TODO +----------- - SSA fabric source/sinks and topo advection terms need to be included +- ... """ import numpy as np +import matplotlib.tri as tri #import code # code.interact(local=locals()) + from ..specfabpy import specfabpy as sf__ +from .. import constants as sfconst from .. import common as sfcom + from firedrake import * class IceFabric: def __init__( - self, mesh, boundaries, L, nu_multiplier=1, nu_realspace=1e-3, modelplane='xz', symframe=-1, ds=None, nvec=None + self, mesh, boundaries, L, *args, + nu_multiplier=1, nu_realspace=1e-3, modelplane='xz', symframe=-1, + Eij_grain=(1,1), alpha=0, n_grain=1, CAFFE_params=(0.1, 10), n_EIE=3, + ds=None, nvec=None, setextra=True, + Cij=sfconst.ice['elastic']['Bennett1968'], rho=sfconst.ice['density'], **kwargs ): ### Check args if modelplane != 'xz': raise ValueError('modelplane "%s" not supported, must be "xz"'%(modelplane)) + + if n_grain != 1: + raise ValueError('only n_grain = 1 (linear viscous) is supported') + + if not(0 <= alpha <= 1): + raise ValueError('alpha should be between 0 and 1') ### Setup + # Firedrake self.mesh, self.boundaries = mesh, boundaries self.ds = ds if ds is not None else Measure('ds', domain=self.mesh, subdomain_data=self.boundaries) self.n = nvec if nvec is not None else FacetNormal(self.mesh) - self.L = int(L) # spectral truncation - self.symframe = symframe - self.modelplane = modelplane - self.nu_realspace = Constant(nu_realspace) - self.nu_multiplier = Constant(nu_multiplier) + self.setextra = setextra # set Eij, E_CAFFE, pfJ, etc. on every state update? - ### Initialize fortran module - + # specfab + self.L = int(L) # spectral truncation self.sf = sf__ self.lm, self.nlm_len = self.sf.init(self.L) self.rnlm_len = self.sf.get_rnlm_len() + self.symframe = symframe + self.modelplane = modelplane + self.nu_realspace = nu_realspace + self.nu_multiplier = nu_multiplier + self.grain_params = (Eij_grain, alpha, n_grain) + self.CAFFE_params = CAFFE_params # (Emin, Emax) + self.n_EIE = n_EIE + self.Lame_grain = self.sf.Cij_to_Lame_tranisotropic(Cij) + self.rho = rho ### Fabric state and dyamics @@ -71,13 +92,15 @@ def __init__( self.R = FunctionSpace(self.mesh, *ele) self.G = TensorFunctionSpace(self.mesh, *ele, shape=(2,2)) self.numdofs = self.R.dim() + # Test, trial, solution container self.p = TrialFunction(self.S) self.w = TestFunction(self.S) -# self.ps = split(self.p) +# self.ps = split(self.p) # split done in _weakform() instead self.ws = split(self.w) self.s = Function(self.S) self.s0 = Function(self.S) + # Dynamical matrices rows, e.g. M_LROT[i,:] (this account for the real-real coefficient interactions, sufficient for xz problems) self.Mrr_LROT = [Function(self.S) for ii in range(self.rnlm_len)] # list of M matrix rows self.Mrr_DDRX_src = [Function(self.S) for ii in range(self.rnlm_len)] @@ -92,11 +115,6 @@ def __init__( self.Gd = TensorFunctionSpace(self.mesh, *ele, shape=(2,2)) self.Vd = VectorFunctionSpace(self.mesh, *ele, dim=3) # for vectors self.numdofs0 = self.Rd.dim() - self.mi = [Function(self.Vd) for _ in range(3)] # (m1,m2,m3) fabric principal directions - self.Eij = [Function(self.Rd) for _ in range(6)] # Eij enhancement tensor - self.lami = [Function(self.Rd) for _ in range(3)] # a2 eigenvalues (lami) - self.E_CAFFE = Function(self.Rd) - self.J = Function(self.Rd) ### Idealized states @@ -117,40 +135,31 @@ def __init__( def initialize(self, s0=None): - s0_ = self.rnlm_iso if s0 is None else s0 - self.s.assign(project(Constant(s0_), self.S)) - self._nlm_nodal() + s0 = self.rnlm_iso if s0 is None else s0 + self.set_state(project(Constant(s0), self.S)) def set_state(self, s): self.s.assign(s) self.s0.assign(s) - + self._setaux() + def set_BCs(self, si, domids): self.bcs = [DirichletBC(self.S, si[ii], did) for ii, did in enumerate(domids)] def set_isotropic_BCs(self, domids): self.set_BCs([Constant(self.rnlm_iso)]*len(domids), domids) - def get_nlm(self, *p, xz2xy=False): - nlm = self.sf.rnlm_to_nlm(self.s(p)+0j, self.nlm_len) - return self.sf.rotate_nlm_xz2xy(nlm) if xz2xy else nlm - - def _nlm_nodal(self): - 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.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]) - - def evolve(self, *args, DDRX_LINEARIZE=2, **kwargs): + def evolve(self, u, *args, DDRX_LINEARIZE=2, **kwargs): """ The fabric solver, called to step fabric field forward in time by amount dt @TODO: this routine and _get_weakform() can surely be optimized by better choice of linear solver, preconditioning, and not assembling the weak form on every solve call """ self.s0.assign(self.s) - F = self._get_weakform(*args, DDRX_LINEARIZE=DDRX_LINEARIZE, **kwargs) + F = self._get_weakform(u, *args, DDRX_LINEARIZE=DDRX_LINEARIZE, **kwargs) FORM = (F==0) if DDRX_LINEARIZE == 0 else (lhs(F)==rhs(F)) solve(FORM, self.s, self.bcs, solver_parameters={'linear_solver':'gmres',}) # non-symmetric system (fastest tested are: gmres, bicgstab, tfqmr) - self._nlm_nodal() + self._setaux(u=u) def _get_weakform(self, u, S, dt, iota=None, Gamma0=None, Lambda0=None, zeta=0, steadystate=False, DDRX_LINEARIZE=2): @@ -188,23 +197,22 @@ def _get_weakform(self, u, S, dt, iota=None, Gamma0=None, Lambda0=None, zeta=0, ### Construct weak form - # Real space advection - F = dot(dot(u, nabla_grad(s)), self.w)*dx + # 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 + F = dot(s_null, self.w)*dx # Time derivative dtinv = Constant(1/dt) if not steadystate: F += dtinv * dot( (s-self.s0), self.w)*dx - # 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 - F += dot(s_null, self.w)*dx + # Real space advection + F += dot(dot(u, nabla_grad(s)), self.w)*dx # Real space stabilization (Laplacian diffusion) - if ENABLE_REG: - F += self.nu_realspace * inner(grad(s), grad(self.w))*dx + F += Constant(self.nu_realspace) * inner(grad(s), grad(self.w))*dx if ENABLE_LROT: F += -sum([ dot(self.Mrr_LROT[ii], s)*self.ws[ii]*dx for ii in self.srrng]) @@ -228,7 +236,8 @@ def _get_weakform(self, u, S, dt, iota=None, Gamma0=None, Lambda0=None, zeta=0, F += F_src - F_snk # Orientation space stabilization (hyper diffusion) - F += -self.nu_multiplier*sum([dot(self.Mrr_REG[ii], s)*self.ws[ii]*dx for ii in self.srrng]) + if ENABLE_REG: + F += -Constant(self.nu_multiplier)*sum([dot(self.Mrr_REG[ii], s)*self.ws[ii]*dx for ii in self.srrng]) return F @@ -250,45 +259,71 @@ def eigenframe(self, *p, **kwargs): mi, lami = sfcom.eigenframe(nlm, symframe=self.symframe, modelplane=self.modelplane) return (mi, lami) - def pfJ(self, *args, **kwargs): + def get_nlm(self, *p, xz2xy=False): + nlm = self.sf.rnlm_to_nlm(self.s(p)+0j, self.nlm_len) + return self.sf.rotate_nlm_xz2xy(nlm) if xz2xy else nlm + + def _setaux(self, u=None): """ - Pole figure J (pfJ) index + 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.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: + self.mi, self.Eij, self.lami = self.get_Eij(ei=()) + self.xi, self.Exij, _ = self.get_Eij(ei=np.eye(3)) + # unpack above eigenframe fields for convenience + self.m1, self.m2, self.m3 = self.mi # rheological symmetry directions + self.E11, self.E22, self.E33, self.E23, self.E31, self.E12 = self.Eij # eigenenhancements + self.Exx, self.Eyy, self.Ezz, self.Eyz, self.Exz, self.Exy = self.Exij # Cartesian enhancements + self.lam1, self.lam2, self.lam3 = self.lami # fabric eigenvalues \lambda_i (not a2 eigenvalues unless symframe=-1) + # additional fields... + self.pfJ = self.get_pfJ() + if u is not None: + self.E_CAFFE = self.get_E_CAFFE(u) + self.E_EIE = self.get_E_EIE(u) - Returns: J + def get_pfJ(self, *args, **kwargs): + """ + Pole figure J (pfJ) index """ - self.J.vector()[:] = sfcom.pfJ(self.nlm, *args, **kwargs)[:] - return self.J + pfJ = Function(self.Rd) + pfJ.vector()[:] = sfcom.pfJ(self.nlm, *args, **kwargs)[:] + return pfJ - def get_E_CAFFE(self, u, Emin=0.1, Emax=10): + def get_E_CAFFE(self, u): """ CAFFE model (Placidi et al., 2010) - - Returns: E_CAFFE """ Df = project(sym(grad(u)), self.Gd).vector()[:] - self.E_CAFFE.vector()[:] = self.sf.E_CAFFE_arr(self.nlm, self.mat3d(Df), Emin, Emax) - return self.E_CAFFE + E_CAFFE = Function(self.Rd) + E_CAFFE.vector()[:] = 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()[:] + E_EIE = Function(self.Rd) +# E_EIE.vector()[:] = self.sf.E_EIE_arr(...) + return E_EIE - def get_Eij(self, Eij_grain, alpha, n_grain, ei=()): + def get_Eij(self, ei=()): """ Bulk enhancement factors wrt ei=(e1,e2,e3) axes for *transversely isotropic* grains. If ei=() then CPO eigenframe is used. - - Returns: (mi, Eij, lami) """ - ### Check args - - if n_grain != 1: - raise ValueError('only n_grain = 1 (linear viscous) is supported') - if not(0 <= alpha <= 1): - raise ValueError('alpha should be between 0 and 1') ### Calculate Eij etc. using specfabpy per node mi, lami = sfcom.eigenframe(self.nlm, symframe=self.symframe, modelplane=self.modelplane) # [node,xyz,i], [node,i] if len(ei) == 0: ei = (mi[:,:,0], mi[:,:,1], mi[:,:,2]) elif len(ei[0].shape) == 1: ei = [np.tile(ei[ii], (self.numdofs0,1)) for ii in range(3)] # ei = (e1[node,3], e2, e3) - Eij = self.sf.Eij_tranisotropic_arr(self.nlm, *ei, Eij_grain, alpha, n_grain) # [node, Voigt vector index 1--6] + Eij = self.sf.Eij_tranisotropic_arr(self.nlm, *ei, *self.grain_params) # [node, Voigt vector index 1--6] # The enhancement factor model depends on effective (homogenized) grain parameters, calibrated against deformation tests. # For CPOs far from the calibration states, negative values *may* occur where Eij should tend to zero if truncation L is not large enough. @@ -296,15 +331,19 @@ def get_Eij(self, Eij_grain, alpha, n_grain, ei=()): ### Populate functions + ei_fs = [Function(self.Vd) for _ in range(3)] # a2 eigenvectors if ei=() + lami_fs = [Function(self.Rd) for _ in range(3)] # a2 eigenvalues (lami) + Eij_fs = [Function(self.Rd) for _ in range(6)] # Eij enhancement tensor + for ii in range(3): # vec 1,2,3 - self.lami[ii].vector()[:] = lami[:,ii] + lami_fs[ii].vector()[:] = lami[:,ii] for jj in range(3): # comp x,y,z - self.mi[ii].sub(jj).vector()[:] = mi[:,jj,ii] + ei_fs[ii].sub(jj).vector()[:] = mi[:,jj,ii] for kk in range(6): - self.Eij[kk].vector()[:] = Eij[:,kk] + Eij_fs[kk].vector()[:] = Eij[:,kk] - return (self.mi, self.Eij, self.lami) + return (ei_fs, Eij_fs, lami_fs) def Gamma0_Lilien(self, u, T, A=4.3e7, Q=3.36e4): """ diff --git a/tests/firedrake/simpleshear-rathmann.py b/tests/firedrake/simpleshear-rathmann.py index 5c0dd7a..dc72316 100644 --- a/tests/firedrake/simpleshear-rathmann.py +++ b/tests/firedrake/simpleshear-rathmann.py @@ -30,14 +30,14 @@ ### Numerics and regularization Nt = 20 # number of time steps to take -L = 6 # spectral truncation *** set L=8 or L=10 unless debugging *** +L = 8 # spectral truncation *** set L=8 or L=10 unless debugging *** -fabric_kwargs = dict(nu_multiplier=1, nu_realspace=1e-4, modelplane='xz') # nu_realspace must be adjusted when changing mesh/resolution (trial-and-error so that solution is stable and smooth) +kwargs_num = dict(nu_multiplier=1, nu_realspace=1e-3, modelplane='xz') # nu_realspace must be adjusted when changing mesh/resolution (trial-and-error so that solution is stable and smooth) ### Fabric dynamics -ENABLE_LROT = 1 # boolean flag -ENABLE_DDRX = 0 # boolean flag +ENABLE_LROT = True +ENABLE_DDRX = False iota = +1 # deck-of-cards behaviour for lattice rotatio #Gamma0 = 1e-1 # uniform DDRX rate factor @@ -45,9 +45,12 @@ ### Viscous anisotropy homogenization parameters -alpha = 0.455 # Taylor--Sachs homogenization weight -Eij_grain = (1, 1e3) # (Ecc, Eca) grain enhancements -n_grain = 1 # grain power-law exponent (only n_grain=1 supported) +alpha = 0.455 # Taylor--Sachs homogenization weight +Eij_grain = (1, 1e3) # (Ecc, Eca) grain enhancements +n_grain = 1 # grain power-law exponent (only n_grain=1 supported) +E_CAFFE = (0.1, 10) # (Emin, Emax) of CAFFE + +kwargs_vaniso = dict(alpha=alpha, Eij_grain=Eij_grain, n_grain=n_grain, E_CAFFE=E_CAFFE) """ Setup firedrake fabric class @@ -56,11 +59,11 @@ ### Mesh and function spaces nx = ny = 16 -mesh = fd.UnitSquareMesh(nx, ny) #, diagonal='crossed') +mesh = fd.UnitSquareMesh(nx, ny, diagonal='right') #, diagonal='crossed') x = fd.SpatialCoordinate(mesh) V = fd.VectorFunctionSpace(mesh, "CG", 1) Q = fd.FunctionSpace(mesh, "CG", 1) # for projecting scalar fabric measures -T = fd.TensorFunctionSpace(mesh, "CG", 1) +T = fd.TensorFunctionSpace(mesh, "CG", 2) ### Velocity field @@ -76,7 +79,7 @@ ### Initialize fabric module boundaries = (1,2,3,4) -fabric = IceFabric(mesh, boundaries, L, **fabric_kwargs) # initializes as isotropic fabric field +fabric = IceFabric(mesh, boundaries, L, **kwargs_num, **kwargs_vaniso) # initializes as isotropic fabric field fabric.set_isotropic_BCs((1,)) # isotropic ice incoming from left-hand boundary, remaining boundaries are free (no fabric fluxes) if not ENABLE_LROT: iota = None @@ -92,13 +95,13 @@ if ENABLE_DDRX: fabric.evolve(u, tau, 30*dt_CFL, iota=iota, Gamma0=Gamma0, steadystate=False, DDRX_LINEARIZE=2) # solving directly for nonlinear steady-state not yet supported, so take a large time step intead to approximate steady-state in plot else: fabric.evolve(u, tau, 1, iota=iota, Gamma0=Gamma0, steadystate=True) -pfJ_steady = fabric.pfJ().copy(deepcopy=True) +pfJ_steady = fabric.get_pfJ().copy(deepcopy=True) """ Time evolution """ -fabric.initialize() # reset to isotropic +fabric.initialize() # reset to isotropic after solving for steady state nn = 0 t = 0.0 @@ -110,9 +113,7 @@ t += dt print("*** Step %i :: dt=%.2e, t=%.2e" % (nn, dt, t)) - fabric.evolve(u, tau, dt, iota=iota, Gamma0=Gamma0) - mi, Eij, lami = fabric.get_Eij(Eij_grain, alpha, n_grain) - E_CAFFE = fabric.get_E_CAFFE(u) + fabric.evolve(u, tau, dt, iota=iota, Gamma0=Gamma0) # automatically updates derived properties (Eij, pfJ, etc.) ### Plot results @@ -136,7 +137,7 @@ ### Plot J index ax = axr1[0] - h = fd.pyplot.tricontourf(fabric.pfJ(), axes=ax, levels=np.arange(1, 3+1e-3, 0.2), extend='max', cmap='YlGnBu') + h = fd.pyplot.tricontourf(fabric.pfJ, axes=ax, levels=np.arange(1, 3+1e-3, 0.2), extend='max', cmap='YlGnBu') cbar = plt.colorbar(h, ax=ax, **kwargs_cb) cbar.ax.set_xlabel(r'$J$ index') h = fd.pyplot.quiver(u, axes=ax, cmap='Reds', width=0.0075) @@ -151,7 +152,7 @@ lvls_E = np.arange(0.6, 2+1e-3, 0.1) divnorm_E = colors.TwoSlopeNorm(vmin=np.amin(lvls_E), vcenter=1, vmax=np.amax(lvls_E)) kwargs_E = dict(levels=lvls_E, norm=divnorm_E, extend='both', cmap='PuOr_r') - h = fd.pyplot.tricontourf(E_CAFFE, axes=ax, **kwargs_E) + h = fd.pyplot.tricontourf(fabric.E_CAFFE, axes=ax, **kwargs_E) cbar = plt.colorbar(h, ax=ax, **kwargs_cb) cbar.ax.set_xlabel(r'$E$ (CAFFE)') @@ -163,7 +164,7 @@ idx = ['11','22','33','23','13','12'] # Voigt ordering for ii, ax in enumerate(axr2[:]): - h = fd.pyplot.tricontourf(Eij[ii], axes=ax, **kwargs_E) + h = fd.pyplot.tricontourf(fabric.Eij[ii], axes=ax, **kwargs_E) cbar = plt.colorbar(h, ax=ax, **kwargs_cb) cbar.ax.set_xlabel(r'$E_{%s}$'%(idx[ii]))