Commit 5106a527 authored by Janosch Preuß's avatar Janosch Preuß
Browse files

cleanup

parent a987bdac
Pipeline #231870 passed with stage
in 10 minutes and 37 seconds
import numpy as np
import matplotlib.pyplot as plt
from colorsys import hsv_to_rgb
from math import pi
from mpmath import pi,sin,cos,exp,log,hankel1,sqrt
Rmin = 3e-1
Rmax = 1e5
def bump(x,z0,D,a):
return exp(log(a)*(D/(x-z0))**2)
def fadeOut(r):
return bump(r,Rmax,Rmax/2,0.5)
def fadeIn(r):
return bump(r,Rmin,8*Rmin,0.5)
def my_color_scheme(z):
if np.isnan(z):
pass
# get polar coordinates, angle in [0,2pi]
phi = np.angle(z)+pi
r = np.absolute(z)
if r <= Rmin:
return hsv_to_rgb(0,0,0),1
if r >= Rmax:
return hsv_to_rgb(0,0,1),1
if np.isnan(z):
h = 0
else:
h = phi/(2*pi)
sat = fadeOut(r)
bri = fadeIn(r)
# Conversion from hsv to rgb colors space
return hsv_to_rgb(h,sat,bri),1
def my_grid_scheme(z):
n_circles = 10
D_circles = 0.5
D_circles = 500
if np.isnan(z):
pass
# get polar coordinates, angle in [0,2pi]
phi = np.angle(z) + pi
r = np.absolute(z)
d_theta = np.min(np.abs( phi/(2*pi) -np.arange(0,n_circles)/n_circles ) )
d_c = np.min(np.abs( r - D_circles*np.arange(0,10*n_circles) ) )
tres_theta = 0.007
#tres_theta = 0.00
tres_c = 40
#if d_theta < tres_theta:
# if d_c > tres_c:
# print("d_theta < tres_theta at but d_c > tres_c z = {0}".format(z))
if d_c < tres_c:
return hsv_to_rgb(0,0,0.0), (tres_c - d_c)/tres_c
elif d_theta < tres_theta:
return hsv_to_rgb(0,0,0.0), (tres_theta - d_theta)/tres_theta
else:
return hsv_to_rgb(0,0,1), 0
def A_over_B(C_a,alpha_a,C_b,alpha_b):
img = np.empty(shape=C_a.shape)
alphas = np.empty(shape=alpha_a.shape)
alphas = alpha_a + alpha_b*(1-alpha_a)
for i in range(C_a.shape[0]):
for j in range(C_a.shape[1]):
img[i][j] = (C_a[i][j]*alpha_a[i,j] + C_b[i][j]*alpha_b[i,j]*(1-alpha_a[i,j] ))/alphas[i,j]
return img,alphas
def plot(f,domain=[-1,1,-1,1],res=(200,200)):
left = domain[0]
right = domain[1]
bottom = domain[2]
top = domain[3]
dx = (right-left)/res[0]
dy = (top-bottom)/res[1]
C_a = np.empty(shape=(res[1],res[0],3))
alpha_a = np.empty(shape=(res[1],res[0]))
C_b = np.empty(shape=(res[1],res[0],3))
alpha_b = np.empty(shape=(res[1],res[0]))
for i in range(res[1]):
y = top - dy*i
for j in range(res[0]):
x = left + dx*j
z = np.complex64(x+1j*y)
fz = f(z)
c_a,a_a = my_color_scheme(fz)
C_a[i][j] = c_a
alpha_a[i,j] = a_a
c_b,a_b = my_grid_scheme(fz)
C_b[i][j] = c_b
alpha_b[i,j] = a_b
#img,alphas = A_over_B(C_a,alpha_a,C_b,alpha_b)
img,alphas = A_over_B(C_b,alpha_b,C_a,alpha_a)
#plt.imshow(C_a,origin="lower",alpha=alpha_a,extent=domain)
#plt.show()
#plt.imshow(C_a,origin="lower",extent=domain)
#plt.show()
plt.imshow(C_a,extent=domain)
plt.show()
plt.imshow(C_b,alpha=alpha_b,extent=domain)
plt.show()
#plt.imshow(C_a,origin="lower",alpha=alpha_a,extent=domain)
plt.imshow(img,alpha=alphas,extent=domain)
plt.xlabel("Re$(\lambda)$")
plt.ylabel("Im$(\lambda)$")
plt.savefig("dtn-VALC-large-window.png",transparent=True)
#plt.imshow(C_a,origin="lower",alpha=alpha_a,extent=domain)
plt.show()
#plt.imshow(img,alpha=alphas,extent=domain)
#plt.axis('off')
#plt.savefig("dtn-VALC-nolabel.png",transparent=True)
#plt.imshow(C_a,origin="lower",alpha=alpha_a,extent=domain)
#plt.show()
from netgen.meshing import *
from netgen.csg import *
from netgen.meshing import Mesh as netmesh
from math import pi
from scipy.sparse import csr_matrix
#from scipy.optimize import minimize,root,least_squares
#import scipy.sparse.linalg
from ngsolve import Mesh as NGSMesh
from compute_pot import local_polynomial_estimator,logder
#np.random.seed(seed=99)
np.random.seed(seed=1)
from ngsolve import *
import ceres_dtn as opt
ngsglobals.msg_level = 0
SetNumThreads(4)
####################################################
# parameters
####################################################
order_ODE = 8 # order of ODE discretization
L_max = 6000 # how many DtN numbers to take
f_hz = 0.007 # frequency
a = 1.0 # coupling radius
L_spacing = 100 # take only DtN numbers with index % L_spacing == 0
spline_order = 1 # order of Spline approx. for coefficients
R_max_ODE = 1.0033 # to how far out the ODE should be solved
Nmax = 6 # maximal dofs for learned IE in radial direction
Nrs = list(range(Nmax))
show_plots = False
nodamping = False
sample_rr = np.linspace(a,R_max_ODE,400)
sample_r = sample_rr.tolist()
l_eval = [10,1000,5000]
eval_lami = [l*(l+1)/a**2 for l in l_eval ]
####################################################
# some auxiliary functions
####################################################
def sqrt_2branch(z):
if z.real < 0 and z.imag <0:
return - sqrt(z)
else:
return sqrt(z)
def P_DoFs(M,test,basis):
return (M[test,:][:,basis])
def Make1DMesh_flex(R_min,R_max,N_elems=50):
mesh_1D = netmesh(dim=1)
pids = []
delta_r = (R_max-R_min)/N_elems
for i in range(N_elems+1):
pids.append (mesh_1D.Add (MeshPoint(Pnt(R_min + delta_r*i, 0, 0))))
n_mesh = len(pids)-1
for i in range(n_mesh):
mesh_1D.Add(Element1D([pids[i],pids[i+1]],index=1))
mesh_1D.Add (Element0D( pids[0], index=1))
mesh_1D.Add (Element0D( pids[n_mesh], index=2))
mesh_1D.SetBCName(0,"left")
mesh_1D.SetBCName(1,"right")
mesh_1D = NGSMesh(mesh_1D)
return mesh_1D
####################################################
# damping model
####################################################
RSun = 6.963e8 # radius of the Sun in m
omega_f = f_hz*2*pi*RSun
gamma0 = 2*pi*4.29*1e-6*RSun
omega0 = 0.003*2*pi*RSun
if f_hz < 0.0053:
gamma = gamma0*(omega_f/omega0)**(5.77)
else:
gamma = gamma0*(0.0053*2*pi*RSun/omega0)**(5.77)
if nodamping:
gamma = 0
#omega = sqrt_2branch(1+2*1j*gamma)*omega_f
omega_squared = omega_f**2+2*1j*omega_f*gamma
####################################################
# load solar models and compute potential
####################################################
rS,cS,rhoS,_,__,___ = np.loadtxt("modelSinput.txt",unpack=True)
rV,vV,rhoV,Tv = np.loadtxt("VALCinput.txt",unpack=True)
# sound speed is not given in VAL-C model
# calculate it from temperature using the ideal gas law
gamma_ideal = 5/3 # adiabatic index
mu_ideal = 1.3*10**(-3) # mean molecular weight
R_gas = 8.31446261815324 # universal gas constant
cV = np.sqrt(gamma_ideal*R_gas*Tv/mu_ideal)
# scale to SI units
cS = (10**-2)*cS
rhoS = (10**3)*rhoS
rhoV = (10**3)*rhoV
# overlap interval
rmin = rV[-1]
rmax = rS[0]
weightS = np.minimum(1, (rmax-rS)/(rmax-rmin))
weightV = np.minimum(1, (rV-rmin)/(rmax-rmin))
r = np.concatenate((rS,rV))
ind = np.argsort(r)
r = r[ind]
c = np.concatenate((cS,cV))[ind]
rho = np.concatenate((rhoS,rhoV))[ind]
weight = np.concatenate((weightS,weightV))[ind]
RSun = 6.963e8 # radius of the Sun in m
c0 = 6.855e5 # sound speed at interface [cm/s]
H = 1.25e7 # density scale height in [cm]
omega_c = 0.0052*2*pi # c0/(2*H) cut-off frequency in Herz
rho_clean = logder(rho,r,weight,1)
c_clean = logder(c,r,weight,1)
i0 = np.min( np.nonzero(r > 0.99)[0] )
f,df,ddf = logder(1/np.sqrt(rho),r,weight,3)
pot_rho = np.sqrt(rho_clean)*(ddf+2*df/r)
pot_c = c_clean
####################################################
# BSpline approximation of c,rho and potential
####################################################
r_beg = r[0]
pot_rho_beg = pot_rho[0]
pot_c_beg = pot_c[0]
r_end = r[-1]
pot_rho_end = pot_rho[-1]
pot_c_end = pot_c[-1]
r = np.append([r_beg for i in range(spline_order)],r)
pot_rho = np.append([pot_rho_beg for i in range(spline_order)],pot_rho)
pot_c = np.append([pot_c_beg for i in range(spline_order)],pot_c)
r = np.append(r,[r_end for i in range(spline_order)])
pot_rho = np.append(pot_rho,[pot_rho_end for i in range(spline_order)])
pot_c = np.append(pot_c,[pot_c_end for i in range(spline_order)])
r = r.tolist()
pot_rho = pot_rho.tolist()
pot_c = pot_c.tolist()
pot_rho_B = BSpline(spline_order,r,pot_rho)(x)
pot_c_B = BSpline(spline_order,r,pot_c)(x)
pot_1D = pot_rho_B - omega_squared/pot_c_B**2
####################################################
# Compute DtN numbers by solving ODEs
####################################################
mesh_1D = Make1DMesh_flex(a,R_max_ODE,100)
sample_r_pts = [mesh_1D(r_pt) for r_pt in sample_r]
fes_DtN = H1(mesh_1D, complex=True, order=order_ODE, dirichlet=[1])
u_DtN,v_DtN = fes_DtN.TnT()
gfu_DtN = GridFunction (fes_DtN)
f_DtN = LinearForm(fes_DtN)
f_DtN.Assemble()
r_DtN = f_DtN.vec.CreateVector()
frees_DtN = [i for i in range(1,fes_DtN.ndof)]
DtN_continuous_eigvals = []
filtered_lam_all = []
DtN_ODE_filtered_all = []
lam_continuous = np.array([l*(l+1)/a**2 for l in np.arange(L_max) ] )
def get_ODE_DtN(lam,DtN_list,save_modes=False):
plot_modes_collect = [sample_rr]
header_modes = "r"
for lami in lam:
a_DtN = BilinearForm (fes_DtN, symmetric=False)
a_DtN += SymbolicBFI((x**2)*grad(u_DtN)*grad(v_DtN))
a_DtN += SymbolicBFI(( lami.item()*a**2 + pot_1D*x**2 )*u_DtN*v_DtN)
a_DtN.Assemble()
gfu_DtN.vec[:] = 0.0
gfu_DtN.Set(1.0, BND)
r_DtN.data = f_DtN.vec - a_DtN.mat * gfu_DtN.vec
gfu_DtN.vec.data += a_DtN.mat.Inverse(freedofs=fes_DtN.FreeDofs()) * r_DtN
if save_modes and lami in eval_lami:
l_idx = eval_lami.index(lami)
eval_gfu = np.array([gfu_DtN(mpt) for mpt in sample_r_pts])
plot_modes_collect.append(eval_gfu.real.copy())
header_modes += " l{0}".format(l_eval[l_idx])
if show_plots:
print("l =", l_eval[l_idx])
plt.plot( sample_rr, eval_gfu.real )
plt.title("Real part of dtn mode l = {0}".format(l_eval[l_idx]))
plt.show()
rows_DtN,cols_DtN,vals_DtN = a_DtN.mat.COO()
A_DtN = csr_matrix((vals_DtN,(rows_DtN,cols_DtN)))
val1 = P_DoFs(A_DtN,[0],frees_DtN).dot(gfu_DtN.vec.FV().NumPy()[frees_DtN])
val2 = P_DoFs(A_DtN,[0],[0]).dot(gfu_DtN.vec.FV().NumPy()[0])
val = -P_DoFs(A_DtN,[0],frees_DtN).dot(gfu_DtN.vec.FV().NumPy()[frees_DtN])-P_DoFs(A_DtN,[0],[0]).dot(gfu_DtN.vec.FV().NumPy()[0])
DtN_list.append(-val.item()/a**2)
# divide by a**2 to get unscaled normal derivative
from cmath import sqrt as csqrt
def lam_to_ell(z):
return -0.5+csqrt( 0.25+a**2*z)
#if z.real > 0:
# return -0.5+csqrt( 0.25+a**2*z)
#else:
# return -0.5-csqrt( 0.25+a**2*z)
def fun(z):
DtN_list = []
get_ODE_DtN([z],DtN_list)
return DtN_list[0]
#domain_dtn = [500*(500+1),3400*(3400+1),-2000000,3000000]
#plot(fun,domain_dtn,res=(700,350))
L_plus = 5000
#L_minus = 6500
L_minus = 5000
#L_minus = 1500
domain_dtn = [-L_minus*(L_minus+1),L_plus*(L_plus+1),-6000000,10000000]
#domain_dtn = [-L_minus*(L_minus+1),L_plus*(L_plus+1),-6000000,19000000]
#domain_dtn = [-3500*(3500+1),5000*(5000+1),-6000000,19000000]
#plot(fun,domain_dtn,res=(500,500))
plot(fun,domain_dtn,res=(500,350))
'''
neg_lls = np.arange(1,L_minus+1)
neg_lambda = -neg_lls*(neg_lls+1)/a**2
pos_lls = np.arange(0,L_plus+1)
pos_lambda = pos_lls*(pos_lls+1)/a**2
lambda_ll = np.append(neg_lambda[::-1],pos_lambda)
dtn_vals = []
for lam_l in lambda_ll:
print("lam_l =" , lam_l)
dtn_vals.append(fun(lam_l))
dtn_vals = np.array(dtn_vals)
plot_collect = [lambda_ll,dtn_vals.real,dtn_vals.imag]
header_str = 'lam real imag'
np.savetxt(fname="dtn-VALC-negL-7.0mHz.dat",
X=np.transpose(plot_collect),
header=header_str,
comments='')
'''
#domain_dtn = [-3500*(3500+1),3400*(3400+1),-2000000,19000000]
#domain_dtn = [domain_dtn[0] ,domain_dtn[1] ,-10000000,19000000]
#domain_dtn = [500*(500+1),3400*(3400+1),-2000000,9000000]
#domain_plot = [-3500*(3500+1),3400*(3400+1),-2000000,19000000]
#plot(fun,domain_dtn,res=(100,100))
####################################################
# Compute poles of dtn
####################################################
'''
from GRPF import Poles_and_Roots
from cmath import sqrt as csqrt
param = {} # for searching poles
param["Tol"] = 1e-9
param["visual"] = 1
param["ItMax"] = 25
param["NodesMax"] = 500000
#param["xrange"] = [-7500*(7500+1),5000*(5000+1)]
#param["xrange"] = [-7500*(7500+1),5000*(5000+1)]
param["xrange"] = [domain_dtn[0],domain_dtn[1]]
param["yrange"] = [domain_dtn[2], domain_dtn[3]]
#param["yrange"] = [-2000000,3000000]
#param["yrange"] = [-2000000,8000000]
param["h"] = 250000
result = Poles_and_Roots(fun,param)
poles = result["poles"]
print("poles of dtn function")
for pole in poles:
print("{0}".format( (pole.real,pole.imag) ) )
print("roots of dtn function")
for root in result["roots"] :
print("{0}".format( (root.real,root.imag) ) )
'''
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment