Commit 79df6f1a authored by Janosch Preuß's avatar Janosch Preuß
Browse files

dtn-VALC plot in intro

parent b53c80e7
Pipeline #231701 passed with stage
in 10 minutes and 58 seconds
......@@ -51,7 +51,7 @@ are submodules of the repository at hand. They can be pulled with
The main tools are:
* [Netgen/NGSove](https://ngsolve.org/) is the FEM library used for implementation in the paper (if you are only interested in the optimization routine for dtn this
* [Netgen/NGSove](https://ngsolve.org/) is the FEM library used for implementation in the thesis (if you are only interested in the optimization routine for dtn this
does not have to be installed). Installation instructions for common operating systems are provided [here](https://ngsolve.org/downloads).
* *ngs_refsol*: Minor extension of *NGSolve* that provides some sample solutions for scattering problems. For installation it should be sufficient to
execute `python3 setup.py install --user` in the top folder of this repository. Alternatively, one may proceed by a manual installation:
......@@ -76,10 +76,10 @@ for the Atmo model. Otherwise, this dependency can be ignored.
We offer two levels of reproduction.
1. The [first](#reproPlots") lets you reproduce the (Latex-based) figures and tables presented in the thesis based on exactly the same data that was used in the thesis.
1. The [first](#reproPlots") lets you reproduce the figures and tables presented in the thesis based on exactly the same data that was used in the thesis.
This data is provided in the folder *data* of this repository. The folder *plots* contains *latex* scripts which can be run to produce the figures and tables.
2. The [second](#reproNumexp) lets you reproduce the data itself. To this end, the scripts contained in the *scripts* folder have to be executed as described below.
If you want to generate the plots based on the data generated by running the scripts simply copy this data into the corresponding *data/...* folder and run the *latex* scripts in the *plots/..*
If you want to generate the plots based on the data generated by running the scripts, simply copy this data into the corresponding *data/...* folder and run the *latex* scripts in the *plots/..*
folder.
# <a name="reproPlots"></a> Reproduction of figures and tables based on already provided data
......@@ -89,7 +89,8 @@ If you want to generate __all__ the Latex-based plots in one go:
cd plot
make figures
Directions to generate only __specific__ figures are given below.
Directions to generate only __specific__ figures (and some non-latex figures) are given below.
The ordering is given according to the chapters of the thesis.
## Chapter 2: DtN maps for time-harmonic waves in separable exterior domains
Change to directory `plots/2_DtN_sep`
......@@ -259,7 +260,7 @@ When the computations have finished the following output data will be contained
* The data for the Atmo model is available in the file *cross-tt-Atmo.dat*.
* The data for the VAL-C model is available in the file *cross-tt-VALC.dat*.
If you would like to reproduce the latex plots of Fig. 6.5 and Fig. 6.6 with the data just computed then
If you would like to reproduce the latex plots of Fig. 6.5 and Fig. 6.6 with the data just computed, then
copy this data from the folder `scripts/6_helio_tbc/obs_helio` into the folder `data/6_helio_tbc/obs_helio` and follow the instructions for Fig. 6.5 and Fig. 6.6 below
for compiling the latex figures.
......@@ -304,16 +305,16 @@ This will produce the following figures and data:
* (__d__): *power-Atmo-N0.png*
* (__e__): *power-Atmo-N1.png*
* (__f__): *power-Atmo-N4.png*
* Fig. 6.9: The data for the first and second horizontal panel is available in files of the form *dtn-Atmo-__RI__-__f__mHz.dat*.
* Fig. 6.9: The data for the first and second horizontal panel is available in files of the form *dtn-Atmo-__RI__- __f__mHz.dat*.
Here, __RI__ is either "real" or "imag", which signifies that this file contains the real or imaginary part of the dtn function, respectively.
The specification __f__ denotes the frequency where [3,5p3,6p5] corresponds to [3.0,5.3,6.5] mHz.
The relative errors shown in the lower panel are contained in the files *dtn-Atmo-relerr-__f__mHz.dat*.
The relative errors shown in the lower panel are contained in the files *dtn-Atmo-relerr- __f__mHz.dat*.
All tables in these files share the same structure.
* The column *l* denotes the index of the eigenvalue.
* The column *ref* is the exact dtn function (solid gray line in figure).
* The columns *nonlocal* and *SHF1a* contain the data for the nonlocal and S-HF-1a condition.
* The columns labelled N__j__ give the results for learned IEs using N=__j__.
* The columns labelled *N__j__* give the results for learned IEs using N=__j__.
The folder `data/6_helio_tbc/tbc_Atmo` contains copies of these data files which you may want to replace with the data just computed. To compile the latex figures,
then change to the directory `plots/6_helio_tbc/tbc_Atmo` and run
......@@ -323,7 +324,7 @@ This will produce the following figures and data:
lualatex dtn-Atmo-6p5mHz.tex
* Fig. 6.10:
* (__a__): The filter is the same as for Fig. 6.13 (___a__) and therefore only output during the reproduction of of the latter figure.
* (__a__): The filter is the same as for Fig. 6.13 (__a__) and therefore only output during the reproduction of of the latter figure.
* (__b__): The plot *time-distance-Atmo-low-filter.png* is directly produced by the script.
* (__c__): The expectation value of the cross covariance is contained in the file *cross-covariance-val-low-theta14.dat*. The first colum *t* denotes the time.
The column *ref* is the reference cross-covariance (solid gray) in the figure. The
......@@ -332,7 +333,7 @@ This will produce the following figures and data:
* (__d__): This is as in (__c__) except that the corresponding filenames are *cross-covariance-val-low-theta28.dat* and *cross-covariance-diff-low-theta28.dat*.
A copy of the produced data is avialable in the folder `data/6_helio_tbc/tbc_Atmo` which you may want to replace with the data just computed. To compile the latex
figures then change to directory `plots/6_helio_tbc/tbc_Atmo` and run
figures, then change to directory `plots/6_helio_tbc/tbc_Atmo` and run
lualatex cross-covariance-low-theta14.tex
lualatex cross-covariance-low-theta14.tex
......@@ -357,12 +358,12 @@ This will produce the following figures and data:
Change to the folder `scripts/6_helio_tbc/tbc_VALC`. Run `python3 power-and-hl-filtered-cross.py "download"`.
The following plots and data products will be produced.
* Fig. 6.12: The plots *power-computed-N_j__-VALC.png* which display the relative errors for N=__j__ with j in [0,1,2,3,4] are directly available after running the script.
* Fig. 6.12: The plots *power-computed-N__j__-VALC.png* which display the relative errors for N=__j__ with j in [0,1,2,3,4] are directly available after running the script.
* Fig 6.13:
* (__a__): The file *phase-filter-low.png* is directly available.
* (__b__): The file *time-distance-low.png* is directly available.
* (__c__): The data for the expectation value of the cross-covariance is stored in the file *cross-covariance-val-low-theta14.dat*. The column *t* is the time, the column *ref* the
reference cross-covariance (solid gray in the figure) and the column *N__j__* are the results with learned IEs using N=__j__ for __j__ in [0,1,2,3,4].
reference cross-covariance (solid gray in the figure) and the columns *N__j__* are the results with learned IEs using N=__j__ for __j__ in [0,1,2,3,4].
The absolute differences shown in the right panel are contained in the file *cross-covariance-diff-low-theta14.dat*.
The folder `data/6_helio_tbc/tbc_VALC` contains a copy of this data which you may want to replace with the results just computed.
To compile the figure, change into the folder `plots/6_helio_tbc/tbc_VALC` and run `lualatex cross-covariance-low-theta.tex`.
......@@ -377,7 +378,7 @@ The following plots and data products will be produced.
### Fig. 6.15
Change to the folder `scripts/6_helio_tbc/tbc_VALC`.
If you executed the instructions to reproduce Fig. 6.12, 6.13, 6.14, then the data is already available and suffices to run `python3 travel-time-accuracy.py`.
If you executed the instructions to reproduce Fig. 6.12, 6.13, 6.14, then the data is already available and it suffices to run `python3 travel-time-accuracy.py`.
Otherwise, run `python3 travel-time-accuracy.py "download"`.
The following data will be produced
......@@ -393,7 +394,7 @@ To compile the figures, change into the folder `plots/6_helio_tbc/tbc_VALC` and
### Fig. 6.16
Change to the folder `scripts/6_helio_tbc/tbc_VALC`.
If you executed the instructions to reproduce Fig. 6.12, 6.13, 6.14 or Fig.14, then the data is already available and suffices to run
If you executed the instructions to reproduce Fig. 6.12, 6.13, 6.14 or Fig.14, then the data is already available and it suffices to run
`python3 double-ridge.py`.
Otherwise, run `python3 double-ridge.py "download"`.
After running the script the figures (__b__) and (__c__) will be available as *double-ridge-VALC-meshed-atmo.png* and *double-ridge-learned-N0.png*.
......@@ -483,7 +484,7 @@ Change to directory `plots/8_sweep_helio/perturb`.
Please switch into the folder `reproduction`.
The scripts for reproducing the experiments are located in subfolders of the folder `scripts`. The subfolders
are labelled according to the corresponding sections of the paper. The scripts are run with the command
`python3 scriptname.py `. While running the scripts the data for reproducing the results of the paper will usually be written to
`python3 scriptname.py `. While running the scripts the data for reproducing the results of the thesis will usually be written to
specific files (located in the folder in which the script is run) which are explained below for each experiment. For comparison, the folder `data` contains the data
which has been used in the thesis.
The ordering below is given according to the chapters of the thesis.
......@@ -494,7 +495,7 @@ Many of the plots in this chapter are obtained as auxiliary data during computat
In these cases we will give specific references to the corresponding reproduction instructions of the subsequent chapters.
### Fig. 2.2
* The real and imaginary parts of dtn are contained in column `zeta` of the files `jump-dtn-kinf16-real.dat` and *jump-dtn-kinf16-imag.dat*, respectively, which are generated while running the script
* The real and imaginary parts of dtn are contained in the column `zeta` of the files *jump-dtn-kinf16-real.dat* and *jump-dtn-kinf16-imag.dat*, respectively, which are generated while running the script
to reproduce [Fig. 4.6](#figJump).
* The poles of dtn, which are shown in the lower panel of the figure as black crosses, are contained in the file *dtn_hom_poles.dat* which is avaialble after running the reproduction instructions for
[Fig 2.3](#figdtnHomExp).
......@@ -536,7 +537,13 @@ Change to the folder `scripts/6_helio_tbc/2_DtN_sep`. Run `python3 solar_coeff.p
Atmo model.
### Fig. 2.8
TODO
Change to the directory `scripts/6_helio_tbc/2_DtN_sep`. Run `python3 dcolor-VALC.py`.
* The values of the dtn function are available in the file *dtn-VALC-ref-0.007Hz.dat*. The column *l* is the index of the eigenvalue.
The columns *zetaReal* and *zetaImag* contains real and imaginary parts of the dtn function, respectively.
* The poles of dtn are available in the file *dtn-VALC-ref-0.007Hz-poles.dat*. The column *poles_real* contains the real part and
the column *poles_imag* the imaginary part.
* The domain coloring plot is output as *dtn-VALC-large-window.png*.
## Chapter 3: Tensor-product discretizations of DtN
Change to directory `scripts/3_DtN_TP`.
......@@ -546,9 +553,9 @@ Change to directory `scripts/3_DtN_TP/3_2_4_SPF/`.
Run `python3 simple_trigpol.py`. Three files are produced.
* (__a__) The file *SPF-trafo_deg6_b1.dat* contains the data for the approximation of the transformed function. The first column *x* is the coordinate.
The second column *ExactReal* is the real part of the exact transformed function at *x* (gray solid line in Fig. 3.3 (a)). The fourth column *AppReal*
The second column *ExactReal* is the real part of the exact transformed function at *x* (gray solid line in Fig. 3.3 (__a__)). The fourth column *AppReal*
is the value of the approximation (blue dashed line).
* (__b__) The file *SPF-original_deg6_b1.dat* contains the data for the approximation of the original function. The columns are labelled as in (a).
* (__b__) The file *SPF-original_deg6_b1.dat* contains the data for the approximation of the original function. The columns are labelled as in (__a__).
* (__c__) The file *SPF-relerr_b1.dat* contains the relative errors. The first column *deg* gives the polynomial degree (N in the figure). The second column *trafoError*
contains the relative error for the transformed function (blue solid line) and the third column *originalError* the relative error for the original function (red solid line) .
......@@ -584,13 +591,13 @@ in Fig. 3.5. The column *lam* is the x-axis so the eigenvalue and the column *co
also contains the approximation by TP-PML (solid orange line) in the column *Nze17*. The file *dtn-U-Burnett-real-omega16.dat* contains the data for U-Burnett (green dashed line) in
the column *N3*. The file *dtn-Astley-Leis-real-omega16.dat* contains the data for Astley-Leis (blue dotted line) in the column *N3*. The file *compTP/dtn-HSIE-real-omega16.dat* contains
the data for the HSIE (red dashdotted) in the column *N3*.
2. panel: The same as in 2. panel except that the specifier *real* in the filenames has to be replaced by *imag*, e.g. has to be used instead of *dtn-TP-PML-real-omega16.dat*.
2. panel: The same as in 1. panel except that the specifier *real* in the filenames has to be replaced by *imag*.
3. panel: The relative error for the TP-PML is in the column *Nze17* of the file *relerrTP-PML-omega16.dat*. The relative error for U-Burnett is in the column *N3* of the file *relerrU-Burnett-omega16.dat*.
The relative error for *Astley-Leis* is in the column *N3* of the file *relerrAstley-Leis-omega16.dat*. The relative error for the HSIE is in the column *N3* of the file *relerrHSIE-omega16.dat*.
4. panel: The exact poles (black crosses) are stored in the file *analytic-poles-omega16.dat*. The column *Re* contains the real parts and the column *Im* the imaginary parts. The columns in the other
files are labelled the same way. The approximate poles for the TP-PML are stored in the file *approx-poles-TP-PML-omega16-Nzes17.dat*. The approximate poles for U-Burnett are stored in
the file *approx-poles-U-Burnett-omega16-N3.dat*. The approximate poles for Astley-Leis are stored in the file *approx-poles-Astley-Leis-omega16-N3.dat*. The approximate poles for the HSIE
are stored in the file *approx-poles-HSIE-omega16-N3.dat.
are stored in the file *approx-poles-HSIE-omega16-N3.dat*.
* (__b__):
1. panel: The same as in (__a__) except that the column *Nze17* for the *TP-PML* needs to be replaced by the column *Nze49* and that the column *N3* for the other methods has to be replaced by *N6*.
2. panel: The same as in (__a__) except that the column *Nze17* for the *TP-PML* needs to be replaced by the column *Nze49* and that the column *N3* for the other methods has to be replaced by *N6*.
......@@ -611,7 +618,7 @@ Change to directory `scripts/4_learnedIE_individual`.
Change to directory `scripts/4_learnedIE_individual/optimization`. Run `python3 learning_comp.py`. Several data files will be generated.
* __Fig. 4.1.__:
* (__a__) The relative errors are stored in files of the form *ansatz-residuals-N__j__.dat* where the integer __j__ denotes the value` of N. Each of these files contains three columns. The first *l*
* (__a__) The relative errors are stored in files of the form *ansatz-residuals-N__j__.dat* where the integer __j__ denotes the value of N. Each of these files contains three columns. The first *l*
describes the index of the eigenvalue as in the figure. The second column *mediumSym* is the relative error for the reduced symmetric ansatz. The last column *full* contains the realtive
errors for the full ansatz (i.e. using dense matrices).
* (__b__) The poles are stored in different files which all have the same structure. Each contains two columns: *poles_real* for the real part of the poles and *poles_imag* for the imaginary part.
......@@ -652,7 +659,7 @@ Change to directory `scripts/4_learnedIE_individual/spherical_geom/point_source`
The following data files are generated:
* For learned IEs based on reduced symmetric ansatz: *pt-source-learned-mediumSym-heu-weights.dat* (blue color in the figure). For learned IEs based on the full ansatz: *pt-source-learned-full-heu-weights.dat* (red color).
* For learned IEs based on the reduced symmetric ansatz: *pt-source-learned-mediumSym-heu-weights.dat* (blue color in the figure). For learned IEs based on the full ansatz: *pt-source-learned-full-heu-weights.dat* (red color).
* For the HSIE: *pt-source-HSIE.dat* (in green color).
* For the TP-PML: *pt-source-TP-PML.dat* (in violet).
* For the adaptive PML: *pt-source-adaptive-PML-close.dat* and *pt-source-adaptive-PML-far.dat* (orange).
......
poles_real poles_imag
1.377503938197673961e+03 1.863660851475341929e+02
2.570556869088773965e+03 1.104691540690125748e+02
3.085935780816633269e+03 8.037360887798470799e+01
......@@ -107,13 +107,15 @@
xmin = 0,
xmax = 6e3,
ymin = -10,
ymax = 200,
ymax = 210,
]
\addplot[green!60!black ,very thick,mark=x, only marks,mark options={scale=3}] coordinates {
(3004.6998574946547,86.00017400788593)
(2582.964868161419,110.37886197860917)
(1493.2648322804205,175.18550306477562)
};
%\addplot[green!60!black ,very thick,mark=x, only marks,mark options={scale=3}] coordinates {
%(3004.6998574946547,86.00017400788593)
%(2582.964868161419,110.37886197860917)
%(1493.2648322804205,175.18550306477562)
%};
\addplot[green!60!black ,very thick,mark=x, only marks,mark options={scale=3}]
table[mark=none,x=poles_real,y=poles_imag] {../../../data/2_DtN_sep/dtn_functions/VALC/dtn-VALC-ref-0.007Hz-poles.dat};
\addplot[mark=none,black,samples=2,very thin,domain=0:6e3] {0.0};
\end{axis}
\begin{axis}[
......@@ -154,7 +156,8 @@
ymax = 20,
]
%\node (mypic) at (64.5,55) {\includegraphics[scale = 0.175]{../VALC-geom.pdf}};
\node (mypic) at (axis cs:4.8,12) {\includegraphics[scale = 0.425]{../../../data/2_DtN_sep/dtn_functions/VALC/dtn-VALC-large-window.png}};
%\node (mypic) at (axis cs:4.8,12) {\includegraphics[scale = 0.425]{../../../data/2_DtN_sep/dtn_functions/VALC/dtn-VALC-large-window.png}};
\node (mypic) at (axis cs:4.8,11.5) {\includegraphics[scale = 0.425]{../../../data/2_DtN_sep/dtn_functions/VALC/dtn-VALC-large-window.png}};
\end{axis}
\end{tikzpicture}
\end{document}
......
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
import sys,os
file_dir = os.path.dirname(os.path.realpath(__file__))
upfolder_path= os.path.abspath(os.path.join( os.path.dirname( __file__ ), '..' ))
sys.path.append(upfolder_path)
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
get_ODE_DtN(lam_continuous,DtN_continuous_eigvals,True)