agd.ExportedCode.Notebooks_NonDiv.NonlinearMonotoneSecond2D

  1# Code automatically exported from notebook NonlinearMonotoneSecond2D.ipynb in directory Notebooks_NonDiv
  2# Do not modify
  3from ... import Selling
  4from ... import LinearParallel as lp
  5from ... import AutomaticDifferentiation as ad
  6from ... import Domain
  7from agd.Plotting import savefig; #savefig.dirName = "Figures/NonlinearMonotoneSecond2D"
  8
  9import numpy as np
 10import matplotlib.pyplot as plt
 11
 12def SchemeNonMonotone(u,α,β,bc,sqrt_relax=1e-6):
 13    # Compute the hessian matrix of u
 14    uxx = bc.Diff2(u,(1,0))
 15    uyy = bc.Diff2(u,(0,1))
 16    uxy = 0.25*(bc.Diff2(u,(1,1)) - bc.Diff2(u,(1,-1)))
 17    
 18    # Compute the eigenvalues
 19    # The relaxation is here to tame the non-differentiability of the square root.
 20    htr = (uxx+uyy)/2. # Half trace
 21    Δ = ((uxx-uyy)/2.)**2 + uxy**2 # Discriminant of characteristic polynomial
 22    sΔ = np.sqrt( np.maximum( Δ, sqrt_relax) )
 23
 24    λ_max = htr+sΔ
 25    λ_min = htr-sΔ
 26    
 27    # Numerical scheme
 28    residue = β - α*λ_max - λ_min
 29    
 30    # Boundary conditions
 31    return np.where(bc.interior,residue,u-bc.grid_values)
 32
 33def SchemeSampling(u,diffs,β,bc):
 34    # Tensor decomposition 
 35    μ,e = Selling.Decomposition(diffs)
 36    
 37    # Numerical scheme 
 38    μ = bc.as_field(μ)
 39    residue = β - (μ*bc.Diff2(u,e)).sum(0).min(0)
 40    
 41    # Boundary conditions
 42    return np.where(bc.interior,residue,u-bc.grid_values)
 43
 44def Diff(α,θ):
 45    e0 = np.array(( np.cos(θ),np.sin(θ)))
 46    e1 = np.array((-np.sin(θ),np.cos(θ)))
 47    if isinstance(α,np.ndarray): 
 48        e0,e1 = (as_field(e,α.shape) for e in (e0,e1))
 49    return α*lp.outer_self(e0) + lp.outer_self(e1)
 50
 51def SchemeSampling_OptInner(u,diffs,bc,oracle=None):
 52    # Select the active tensors, if they are known
 53    if not(oracle is None):
 54        diffs = np.take_along_axis(diffs, np.broadcast_to(oracle,diffs.shape[:2]+(1,)+oracle.shape),axis=2)
 55    
 56    print("Has AD information :", ad.is_ad(u), ". Number active tensors per point :", diffs.shape[2])
 57    
 58    # Tensor decomposition 
 59    coefs,offsets = Selling.Decomposition(diffs)
 60    
 61    # Return the minimal value, and the minimizing index
 62    return ad.min_argmin( lp.dot_VV(coefs,bc.Diff2(u,offsets)), axis=0)
 63
 64def SchemeSampling_Opt(u,diffs,β,bc):
 65    # Evaluate the operator using the envelope theorem
 66    result,_ = ad.apply(SchemeSampling_OptInner, u,bc.as_field(diffs),bc, envelope=True)
 67        
 68    # Boundary conditions
 69    return np.where(bc.interior, β-result, u-bc.grid_values)
 70
 71def MakeD(α):
 72    return np.moveaxis(0.5*np.array([
 73        (α+1)*np.array([[1,0],[0,1]]),
 74        (α-1)*np.array([[1,0],[0,-1]]),
 75        (α-1)*np.array([[0,1],[1,0]])
 76    ]), 0,-1)
 77
 78def NextAngleAndSuperbase(θ,sb,D):
 79    pairs = np.stack([(1,2), (2,0), (0,1)],axis=1)
 80    scals = lp.dot_VAV(np.expand_dims(sb[:,pairs[0]],axis=1), 
 81                       np.expand_dims(D,axis=-1), np.expand_dims(sb[:,pairs[1]],axis=1))
 82    ϕ = np.arctan2(scals[2],scals[1])
 83    cst = -scals[0]/np.sqrt(scals[1]**2+scals[2]**2)
 84    θ_max = np.pi*np.ones(3)
 85    mask = cst<1
 86    θ_max[mask] = (ϕ[mask]-np.arccos(cst[mask]))/2
 87    θ_max[θ_max<=0] += np.pi
 88    θ_max[θ_max<=θ] = np.pi
 89    k = np.argmin(θ_max)
 90    i,j = (k+1)%3,(k+2)%3
 91    return (θ_max[k],np.stack([sb[:,i],-sb[:,j],sb[:,j]-sb[:,i]],axis=1))
 92
 93def AnglesAndSuperbases(D,maxiter=200):
 94    sb = Selling.CanonicalSuperbase(np.eye(2)).astype(int)
 95    θs=[]
 96    superbases=[]
 97    θ=0
 98    for i in range(maxiter):
 99        θs.append(θ)
100        if(θ>=np.pi): break
101        superbases.append(sb)
102        θ,sb = NextAngleAndSuperbase(θ,sb,D)
103    return np.array(θs), np.stack(superbases,axis=2)
104
105def MinimizeTrace(u,α,bc,sqrt_relax=1e-16):
106    # Compute the tensor decompositions
107    D=MakeD(α)
108    θ,sb = AnglesAndSuperbases(D)
109    θ = np.array([θ[:-1],θ[1:]])
110    
111    # Compute the second order differences in the direction orthogonal to the superbase
112    sb_rotated = np.array([-sb[1],sb[0]])
113    d2u = bc.Diff2(u,sb_rotated)
114    d2u[...,bc.not_interior]=0. # Placeholder values to silent NaNs
115    
116    # Compute the coefficients of the tensor decompositions
117    sb1,sb2 = np.roll(sb,1,axis=1), np.roll(sb,2,axis=1)
118    sb1,sb2 = (e.reshape( (2,3,1)+sb.shape[2:]) for e in (sb1,sb2))
119    D = D.reshape((2,2,1,3,1)+D.shape[3:])
120    # Axes of D are space,space,index of superbase element, index of D, index of superbase, and possibly shape of u
121    scals = lp.dot_VAV(sb1,D,sb2)
122
123    # Compute the coefficients of the trigonometric polynomial
124    scals,θ = (bc.as_field(e) for e in (scals,θ))
125    coefs = -lp.dot_VV(scals, np.expand_dims(d2u,axis=1))
126    
127    # Optimality condition for the trigonometric polynomial in the interior
128    value = coefs[0] - np.sqrt(np.maximum(coefs[1]**2+coefs[2]**2,sqrt_relax))
129    coefs_ = ad.remove_ad(coefs) # removed AD information
130    angle = np.arctan2(-coefs_[2],-coefs_[1])/2.
131    angle[angle<0]+=np.pi
132    
133    # Boundary conditions for the trigonometric polynomial minimization
134    mask = np.logical_not(np.logical_and(θ[0]<=angle,angle<=θ[1]))
135    t,c = θ[:,mask],coefs[:,mask]
136    value[mask],amin_t = ad.min_argmin(c[0]+c[1]*np.cos(2*t)+c[2]*np.sin(2*t),axis=0)
137        
138    # Minimize over superbases
139    value,amin_sb = ad.min_argmin(value,axis=0)
140    
141    # Record the optimal angles for future use
142    angle[mask]=np.take_along_axis(t,np.expand_dims(amin_t,axis=0),axis=0).squeeze(axis=0) # Min over bc
143    angle = np.take_along_axis(angle,np.expand_dims(amin_sb,axis=0),axis=0) # Min over superbases
144
145    return value,angle
146
147def SchemeConsistent(u,α,β,bc):
148    value,_ = MinimizeTrace(u,α,bc)
149    residue = β - value
150    return np.where(bc.interior,residue,u-bc.grid_values)
151
152def MinimizeTrace_Opt(u,α,bc,oracle=None):
153    if oracle is None:  return MinimizeTrace(u,α,bc)
154    
155    # The oracle contains the optimal angles
156    diffs=Diff(α,oracle.squeeze(axis=0))
157    coefs,sb = Selling.Decomposition(diffs)
158    value = lp.dot_VV(coefs,bc.Diff2(u,sb))
159    return value,oracle
160    
161
162def SchemeConsistent_Opt(u,α,β,bc):
163    value,_ = ad.apply(MinimizeTrace_Opt,u,α,bc,envelope=True)
164    residue = β - value
165    return np.where(bc.interior,residue,u-bc.grid_values)
166
167def Pucci_ad(u,α,x):
168    """
169    Computes alpha*lambda_max(D^2 u) + lambda_min(D^2 u), 
170    at the given set of points, by automatic differentiation.
171    """
172    x_ad = ad.Dense2.identity(constant=x,shape_free=(2,))
173    hessian = u(x_ad).hessian()
174    
175    Δ = ((hessian[0,0]-hessian[1,1])/2.)**2 + hessian[0,1]**2
176    sΔ = np.sqrt(Δ)
177    mean = (hessian[0,0]+hessian[1,1])/2.
178    λMin,λMax = mean-sΔ,mean+sΔ
179    
180    return λMin+α*λMax