Source code for nse_extlin_utils

import numpy as np
import scipy.linalg as spla
import dolfin_navier_scipy.stokes_navier_utils as snu
import dolfin_navier_scipy.data_output_utils as dou


[docs]def get_get_cur_extlin(vinf=None, V=None, diribcs=None, invinds=None, reducedmodel=False, tl=None, tr=None, picrdvsnwtn=0., amat=None, akmat=None, **kwargs): ''' returns a function to compute the current extended linearization Parameters --- vinf: array the set point use_ric_ini : string, optional path to a stabilizing initial guess ''' convc_mat, vinablavi, _ = \ snu.get_v_conv_conts(prev_v=vinf, invinds=invinds, retparts=False, V=V, diribcs=diribcs) def get_cur_extlin(vcur=None, vdelta=None): ''' compute the current extended linearization dot vd + A(vd)vd = rhs for the NSE momentum eqn (shifted by the setpoint vinf) which basically reads A(vd)(.) = amat*(.) + (.*nabla)*vinf + ([vinf+vd]*nabla)*(.) or A(vd)(.) = amat*(.) + (.*nabla)*[vinf+vd] + (vinf*nabla)*(.) or anything in between Parameters --- vdelta: array the difference between the (target) set point and the current state vcur: nparray the current state picrdvsnwtn: float, optional blends between the picard (0.) and antipicard (1.) part of the extended linearization ''' if vdelta is None: vdelta = vcur - vinf fullnv = V.dim() vdpbc = np.zeros((fullnv, 1)) # vdelta has zero bcs vdpbc[invinds] = vdelta delta_convmats, vdnablavd, _ = \ snu.get_v_conv_conts(prev_v=vdpbc, invinds=invinds, retparts=True, zerodiribcs=True, V=V, diribcs=None) # testing to see the effect of the zero dirichlets... # tstdelta_convmats, tstvdnablavd, tstrhs = \ # snu.get_v_conv_conts(prev_v=vdelta, invinds=invinds, # retparts=True, # V=V, diribcs=diribcs) # delta_convmats[0] -- Picard part # delta_convmats[1] -- Anti Picard (Newton minus Picard) # one may think of convex combining both parts if picrdvsnwtn < 0 or picrdvsnwtn > 1: raise UserWarning('interesting parameter value -- good luck!!!') curfmat = amat + convc_mat + \ picrdvsnwtn*delta_convmats[0] + (1-picrdvsnwtn)*delta_convmats[1] return curfmat return get_cur_extlin
def _pinsvd(jmat=None, Minv=None): Np, Nv = jmat.shape minvjt = Minv(jmat.T).todense() S = jmat*minvjt Sinv = spla.inv(S) Pi = np.eye(Nv) - minvjt.dot(Sinv*jmat) umat, svec, vmat = spla.svd(Pi) uk = umat[:, :Nv-Np] vk = vmat.T[:, :Nv-Np] thl = uk*svec[:Nv-Np] thr = vk return thl, thr def decomp_leray(mmat, jmat, dataprfx='', testit=False): from spacetime_galerkin_pod.ldfnp_ext_cholmod import SparseFactorMassmat facmy = SparseFactorMassmat(mmat) Minv = facmy.solve_M pinsvddict = dict(jmat=jmat, Minv=Minv) thrthlstrl = [dataprfx+'_pithl', dataprfx+'_pithr'] thl, thr = dou.load_or_comp(arraytype='dense', filestr=thrthlstrl, numthings=2, comprtn=_pinsvd, comprtnargs=pinsvddict) minvthr = Minv(thr) if testit: Np, Nv = jmat.shape print(np.allclose(np.eye(Nv-Np), thr.T.dot(thl))) minvjt = Minv(jmat.T).todense() S = jmat*minvjt Sinv = spla.inv(S) Pi = np.eye(Nv) - minvjt.dot(Sinv*jmat) print(np.allclose(Pi, thl.dot(thr.T))) return thl, thr, minvthr