agd.ExportedCode.Notebooks_FMM.HighAccuracy 

 1# Code automatically exported from notebook HighAccuracy.ipynb in directory Notebooks_FMM
 2# Do not modify
 3from ... import Eikonal
 4from ... import Metrics
 5from agd.Metrics.Seismic import Hooke
 6from ... import FiniteDifferences as fd
 7from ... import LinearParallel as lp
 8from ... import AutomaticDifferentiation as ad
 9from agd.Plotting import savefig; #savefig.dirName = 'Figures/HighAccuracy'
10
11import numpy as np
12import matplotlib.pyplot as plt
13
14def PoincareCost(q):
15    """
16    Cost function defining the Poincare half plane model of the hyperbolic plane.
17    """
18    return 1/q[1]
19
20def PoincareDistance(p,q):
21    """
22    Distance between two points of the half plane model of the hyperbolic plane.
23    """
24    a = p[0]-q[0]
25    b = p[1]-q[1]
26    c = p[1]+q[1]
27    d = np.sqrt(a**2+b**2)
28    e = np.sqrt(a**2+c**2)
29    return np.log((e+d)/(e-d))
30
31diagCoef = (0.5**2,1) # Diagonal coefficients of M
32
33def diff(x,y,α=0.5): return ad.array([x,y+α*np.sin(np.pi*x)]) # Diffeomorphism f
34
35def RiemannMetric(diag,diff,x,y): 
36    X_ad = ad.Dense.identity(constant=(x,y),shape_free=(2,))
37    Jac = np.moveaxis(diff(*X_ad).gradient(),0,1)
38    return Metrics.Riemann.from_diagonal(diag).inv_transform(Jac)
39
40def RiemannExact(diag,diff,x,y):
41    a,b = diag
42    fx,fy = diff(x,y)
43    return np.sqrt(a*fx**2+b*fy**2)
44
45M=((1.25,0.5),(0.5,2.))
46
47def v(x,y,γ): return γ*np.sin(np.pi*x)*np.sin(np.pi*y)/np.pi
48
49def RanderMetric(x,y,γ=0.8):
50    X_ad = ad.Dense.identity(constant=(x,y),shape_free=(2,))
51    omega = v(*X_ad,γ).gradient()
52    return Metrics.Rander(M,omega)
53
54def RanderSolution(x,y,γ=0.8):
55    return Metrics.Riemann(M).norm((x,y)) + v(x,y,γ)
56
57def ConformalMap(x):
58    """
59    Implements the mapping x -> (1/2) * x^2, where x is seen as a complex variable.
60    """
61    return ad.array([0.5*(x[0]**2-x[1]**2), x[0]*x[1]])
62
63def ConformalApply(norm,f,x,decomp=True):
64    """
65    Applies a conformal change of coordinates to a norm.
66    decomp : decompose the Jacobian into a scaling and rotation.
67    """
68    x_ad = ad.Dense.identity(constant=x,shape_free=(2,))
69    Jac = np.moveaxis(f(x_ad).gradient(),0,1)
70    if not decomp: return norm.inv_transform(Jac)
71    
72    # Assuming Jac is a homothety-rotation
73    α = np.power(lp.det(Jac), 1/norm.vdim)
74    R = Jac/α
75    return norm.with_cost(α).rotate(lp.transpose(R))
76    
77
78def MappedNormValues(norm,f,x,seed):
79    seed = fd.as_field(seed,x.shape[1:])
80    return norm.norm(f(x)-f(seed))
def PoincareCost(q): 
15def PoincareCost(q):
16    """
17    Cost function defining the Poincare half plane model of the hyperbolic plane.
18    """
19    return 1/q[1]

Cost function defining the Poincare half plane model of the hyperbolic plane.

def PoincareDistance(p, q): 
21def PoincareDistance(p,q):
22    """
23    Distance between two points of the half plane model of the hyperbolic plane.
24    """
25    a = p[0]-q[0]
26    b = p[1]-q[1]
27    c = p[1]+q[1]
28    d = np.sqrt(a**2+b**2)
29    e = np.sqrt(a**2+c**2)
30    return np.log((e+d)/(e-d))

Distance between two points of the half plane model of the hyperbolic plane.

diagCoef = (0.25, 1)
def diff(x, y, α=0.5): 
34def diff(x,y,α=0.5): return ad.array([x,y+α*np.sin(np.pi*x)]) # Diffeomorphism f
def RiemannMetric(diag, diff, x, y): 
36def RiemannMetric(diag,diff,x,y): 
37    X_ad = ad.Dense.identity(constant=(x,y),shape_free=(2,))
38    Jac = np.moveaxis(diff(*X_ad).gradient(),0,1)
39    return Metrics.Riemann.from_diagonal(diag).inv_transform(Jac)
def RiemannExact(diag, diff, x, y): 
41def RiemannExact(diag,diff,x,y):
42    a,b = diag
43    fx,fy = diff(x,y)
44    return np.sqrt(a*fx**2+b*fy**2)
M = ((1.25, 0.5), (0.5, 2.0))
def v(x, y, γ): 
48def v(x,y,γ): return γ*np.sin(np.pi*x)*np.sin(np.pi*y)/np.pi
def RanderMetric(x, y, γ=0.8): 
50def RanderMetric(x,y,γ=0.8):
51    X_ad = ad.Dense.identity(constant=(x,y),shape_free=(2,))
52    omega = v(*X_ad,γ).gradient()
53    return Metrics.Rander(M,omega)
def RanderSolution(x, y, γ=0.8): 
55def RanderSolution(x,y,γ=0.8):
56    return Metrics.Riemann(M).norm((x,y)) + v(x,y,γ)
def ConformalMap(x): 
58def ConformalMap(x):
59    """
60    Implements the mapping x -> (1/2) * x^2, where x is seen as a complex variable.
61    """
62    return ad.array([0.5*(x[0]**2-x[1]**2), x[0]*x[1]])

Implements the mapping x -> (1/2) * x^2, where x is seen as a complex variable.

def ConformalApply(norm, f, x, decomp=True): 
64def ConformalApply(norm,f,x,decomp=True):
65    """
66    Applies a conformal change of coordinates to a norm.
67    decomp : decompose the Jacobian into a scaling and rotation.
68    """
69    x_ad = ad.Dense.identity(constant=x,shape_free=(2,))
70    Jac = np.moveaxis(f(x_ad).gradient(),0,1)
71    if not decomp: return norm.inv_transform(Jac)
72    
73    # Assuming Jac is a homothety-rotation
74    α = np.power(lp.det(Jac), 1/norm.vdim)
75    R = Jac/α
76    return norm.with_cost(α).rotate(lp.transpose(R))

Applies a conformal change of coordinates to a norm. decomp : decompose the Jacobian into a scaling and rotation.

def MappedNormValues(norm, f, x, seed): 
79def MappedNormValues(norm,f,x,seed):
80    seed = fd.as_field(seed,x.shape[1:])
81    return norm.norm(f(x)-f(seed))