n-claes / legolas Goto Github PK
View Code? Open in Web Editor NEWLegolas: a modern tool for MHD spectroscopy
Home Page: https://legolas.science
License: GNU General Public License v3.0
Legolas: a modern tool for MHD spectroscopy
Home Page: https://legolas.science
License: GNU General Public License v3.0
On the using pylbo page on the website, the link "installing pylbo" below the table of contents is broken.
When enabling both the continuum and eigenfunction plotting features these interfere with each other. Clicking on the legend also selects the nearest spectrum point, which should not happen.
Possible solutions:
Low priority for now
Add boundary conditions when calculating array derivatives in get_array_derivative
subroutine. As of right now values of neighbouring cells are simply copied to the edges. Either find a way to consistently hardcode this somehow, or ask for user-defined input.
After a recent push to one of the development branches the docs are suddenly failing:
Traceback (most recent call last):
[8](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:9)
File "/opt/hostedtoolcache/Python/3.9.13/x64/bin/ford", line 8, in <module>
[9](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:10)
sys.exit(run())
[10](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:11)
File "/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/ford/__init__.py", line 633, in run
[11](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:12)
proj_data, proj_docs, md = initialize()
[12](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:13)
File "/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/ford/__init__.py", line 236, in initialize
[13](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:14)
return parse_arguments(vars(args), proj_docs, directory)
[14](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:15)
File "/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/ford/__init__.py", line 406, in parse_arguments
[15](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:16)
md = markdown.Markdown(
[16](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:17)
File "/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/markdown/core.py", line 96, in __init__
[17](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:18)
self.registerExtensions(extensions=kwargs.get('extensions', []),
[18](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:19)
File "/opt/hostedtoolcache/Python/3.9.13/x64/lib/python3.9/site-packages/markdown/core.py", line 125, in registerExtensions
[19](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:20)
ext.extendMarkdown(self)
[20](https://github.com/n-claes/legolas/runs/7408674995?check_suite_focus=true#step:5:21)
TypeError: extendMarkdown() missing 1 required positional argument: 'md_globals'
This seems to be an upstream issue with Ford...
Performing multiruns for the interchange modes yields figures that seem 'flipped' with respect to the expected values. Perhaps this is also due to the exponential behaviour, but more likely it's a bug in setting up the equilibrium itself.
Since the removal of the grid(1)
inconsistency in cylindrical coordinates in #13
(commit d978dba), the first value of multiple eigenfunctions (that is, the value corresponding to the first gridpoint) is NaN
, because we rescale the values by dividing by r
(which is 0).
This is probably the one and only reason why grid(1)
was set to 1e-5
in the original LEDA code. I think it should be possible to circumvent this by restructuring how the eigenfunctions are composed from the eigenvectors, but it's best to wait on #16 before doing this (so it can be integrated directly).
There seems to be a small discrepancy in calculating the full resistivity and thermal conduction.
pi
instead of 4sqrt(2pi)/3
, which differs by approximately 0.2. Also, log(coulomb_log)
should actually be just coulomb_log
.This has gone unnoticed up to now, since we manually forced the resistivity to a specific value in all currently implemented equilibria. Also the perpendicular thermal conduction coefficient was mostly set to 0 up to now, or is just so small that it has almost no effect.
Nevertheless, implementation should be correct and we should fix this.
At the moment, pylbo handles datasets as follows:
A typical use case handles a few eigenfunctions at most, such that loading all eigenfunctions into memory is really inefficient and a waste of resources.
My current idea to improve this:
I have to test this in more detail, but I'm confident this will have a major impact on memory management and loading times.
In one specific case the setup script setuplegolas.py
ignores the user submodule. When "no" is answered to the question if the default template should be copied over, and smod_user_defined
is added afterwards, the build process ignores the user file and defaults to the one in the main repository.
Minimal example for reproduction
setuplegolas.py
>> No user-defined submodule found, copy default template? [y/n] n
>> Submodule not copied, using default in Legolas src directory.
Add the file, change the equilibrium to "user_defined"
.
Add the cooling curve from Rosner, Tucker & Vaiana (1978) to the code. Don't use interpolation for this one, but rely on the analytical expression.
Multiruns for the coronal/photospheric flux tubes from Bernard Robert's book show a different scaling compared to the the theoretical figures from the book.
Coronal flux tube: theoretical values start at 0.90, legolas values start at 0.995
Photospheric flux tube: theoretical values start a little under 0.90, legolas values start at 0.95
Strange behaviour occurs when the Alfvén or slow continua are plotted versus the grid. There seems to be a discontinuity in some cases (eg. for the resistive tearing modes) which should not occur, these should be smooth curves.
Not sure yet what is going wrong here...
The search function (magnifying glass icon) on the website returns no search results.
Pylbo's set_scaling
method is inconsistent when multiple datasets are added to the same figure. When more than one scaling factor is applied to the datasets the last one always takes precedence and overrides.
Minimal example for reproduction
import matplotlib.pyplot as plt
import pylbo
fig, ax = plt.subplots(1)
series1 = pylbo.load_series(files_1)
series2 = pylbo.load_series(files_2)
p1 = pylbo.plot_spectrum_multi(series1, xdata="k3", custom_figure=(fig, ax))
p1.set_y_scaling(0.5)
p2 = pylbo.plot_spectrum_multi(series2, xdata="k3", custom_figure=(fig, ax))
p2.set_y_scaling(0.1)
Actual result
The scaling factor of the first dataset is overridden by 0.1.
Expected result
Separate scaling factors for both datasets.
Quadratic and cubic functions seem incorrect, something strange is going on with rj_lo
, rj
and rj_hi
. Seems that we only have to use rj_lo
and rj_hi
but this needs to be figured out. Implementation is kept as-is for the moment.
Make sure that all submodules have extensive docs + regression tests, and that Pylbo has unit tests.
Current progress:
When using plot_merged_spectrum
in combination with add_eigenfunctions
to combine datfiles with only 1 eigenvalue, the eigenfunctions cannot be displayed due to an index error.
Minimal example for reproduction
The passed argument is a folder path to a folder containing datfiles with a single eigenvalue and its associated eigenfunctions.
from os import listdir
import pylbo
import argparse as ap
parser = ap.ArgumentParser()
parser.add_argument('-i', required=True, action='store', dest='i', help='followed by file path')
parser.add_argument('-ef', action='store_true', default=False, help='show eigenfunctions')
args = parser.parse_args()
files = listdir(args.i)
for ix in range(len(files)):
files[ix] = args.i + files[ix]
files.sort()
series = pylbo.load_series(files)
p = pylbo.plot_merged_spectrum(series)
if args.ef:
p.add_eigenfunctions()
p.show()
Actual result
Traceback (most recent call last):
File "/home/jordi/anaconda3/envs/legolas/lib/python3.8/site-packages/matplotlib/cbook/__init__.py", line 287, in process
func(*args, **kwargs)
File "/home/jordi/codes/legolas/post_processing/pylbo/visualisation/eigenfunction_interface.py", line 174, in on_point_pick
self.on_left_click(event)
File "/home/jordi/codes/legolas/post_processing/pylbo/visualisation/eigenfunction_interface.py", line 219, in on_left_click
idx, xdata, ydata = self._get_clicked_point_data(event)
File "/home/jordi/codes/legolas/post_processing/pylbo/visualisation/eigenfunction_interface.py", line 299, in _get_clicked_point_data
return idx, xdata[idx], ydata[idx]
File "/home/jordi/anaconda3/envs/legolas/lib/python3.8/site-packages/numpy/ma/core.py", line 3222, in __getitem__
dout = self.data[indx]
IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed
Expected result
Display of the eigenfunctions of the requested mode.
Version info
In the new Arnoldi rework #97 we solve a linear system AX = B
using LAPACK's zgbsv routine. This works quite well for now, but this recreates the LU-decomposition at every iteration. Performance can be (greatly?) improved by calculating the LU-decomposition before the iteration and using it at every step.
This will change the solve_linear_system_complex_banded
routine into two calls:
As PR #97 introduced a lot of changes in the solver routines already, we should follow this order:
develop
develop
and then mergedevelop
and create the PR with this optimisationThat way the new changes can directly benefit from the updated regression framework.
We keep extending the submodules containing new equilibria, so it's maybe a good idea to further divide them in subfolders (resistivity, flow, conduction, etc). That way they remain sorted and we don't loose overview of what's what.
Perhaps do this after #25 and when the project board tests & checks
has less 'To do' things.
Extend submodules with additional equilibria containing radiative cooling and conduction. Start with results from Keppens & Van Der Linden, Solar Physics 144, 1993. This paper uses the cooling curve by Rosner et al., so #24 should be done first.
The constant_current_m-2
multirun shows almost the same result as expected, however, there are two additional branches present in the bottom-right corner which should not appear for this case.
Not sure what is going wrong here...
2 issues with this:
1. Arpack solvers
Right now an Arpack-based solver only solves for a part of the total amount of eigenvalues present in the matrix, but saves the whole eigenvalue and eigenfunction array in memory + to the datfile.
Especially for large resolutions this results in a lot of wasted data usage: say we solve for 100/1000 eigenvalues with eigenfunctions, this means that 900 eigenvalues are simply set to 0 and saved to disk, and _all_their eigenfunctions are automatically set to zero -- hence an (ef_gridpts x 900) matrix of zeros wasted both in terms of memory and storage!!
2. Left eigenvectors
We also allocate the left eigenvectors, but these are never used (and hence also take up memory).
Possible solution
Only allocate arrays that are actually used. In principle this shouldn't be too difficult to implement, these are things that should be taken into account:
All other stuff should not be affected, provided those subroutines use inferred arraysizes (which should be the case almost everywhere).
For gcc <= 8 the current develop
branch throws a segmentation fault at runtime (compilation is fine). From gcc 9 onwards everything runs fine.
Minimal example for reproduction
The problematic piece of code is a polymorphic assignment like this
class(*), allocatable :: element
element = (1.0d0, 1.0d0)
! or
allocate(element, source=(1.0d0, 1.0d0))
Both these lines throw a segmentation fault. May be similar to this gcc bug.
I've traced the problem down to the polymorphic functions get_diagonal_factor
and get_quadblock_element
in the essential boundaries and natural boundaries, respectively.
The implementation of the natural boundary conditions is not correct. For most of the "wall" boundary conditions this does not matter, since the relevant variables are zeroed out at the boundary anyway. However, in some cases (eg. additional terms in the energy equation if resistivity is taken into account) these extra terms should be added.
As of right now there are lots of files saved each run: the eigenvalues, grid, equilibrium + names, eigenfunctions etc, resulting in about 10-15 files per run if everything is written.
This is quite cumbersome also maintainance-wise, so my idea is to store everything (except the matrices) into one single .dat file, where we write a header first which stores the offsets for every 'type' of data (similarly as to how MPI-AMRVAC writes snapshots).
Everything can then be read in using Python's struct.unpack()
using those offsets. We'll have to write a specifically-tailored file reader, but in the long run this would be much more beneficial. A format like this is also easily extendible.
Bug report: after running inverse vector-iteration, the result is one eigenmode. When selecting this one eigenmode, an error is thrown because eigenfunc_interface._get_clicked_point_data tries to apply an index to something that is no list.
Minimal example for reproduction
#Parfile:
&gridlist
gridpoints = 100
/
&equilibriumlist
equilibrium_type = "MRI_accretion"
use_defaults = .false.
/
&savelist
write_eigenfunctions = .true.
write_derived_eigenfunctions = .false.
show_results = .true.
basename_datfile = "sari_GK2022_quasi_continuum_test"
/
&physicslist
flow = .true.
external_gravity = .true.
incompressible = .false.
/
¶mlist
k2 = 1.0d0
k3 = 50.0d0
delta = 0.28862653d0 ! for most unstable mode
beta = 10.0d0
tau = 10.0d0 !fixes mu1 parameter
nu = 0.045d0 !fixes epsilon parameter
/
&solvelist
solver = "inverse-iteration"
sigma = (0.59639333d0, 0.17629975d0)
maxiter = 20000
tolerance = 1.d-7
/
Actual result
File "/users/cpa/nicolasb/codes/legolas/post_processing/pylbo/visualisation/eigenfunctions/eigfunc_interface.py", line 262, in on_left_click
idx, xdata, ydata = self._get_clicked_point_data(event)
File "/users/cpa/nicolasb/codes/legolas/post_processing/pylbo/visualisation/eigenfunctions/eigfunc_interface.py", line 345, in _get_clicked_point_data
return idx, xdata[idx], ydata[idx]
File "/users/cpa/nicolasb/Documents/no_backup/miniconda3/envs/legolasenv/lib/python3.9/site-packages/numpy/ma/core.py", line 3224, in __getitem__
dout = self.data[indx]
IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed
Expected result
The eigenvalue is selected and everyone lives happily ever after.
When selecting xdata, ydata from the artist in _get_cicked_point_data, the result is probably a float instead of a list now if the spectrum consists of only one point.
Version info
With future applications in mind, look into a parallelisation of the BLAS and LAPACK libraries. In essence this is simply a call to eg. SCALAPACK, so we have to figure out how to set that up during the make
process. Another option is using the Intel MKL parallel libraries for this.
Low priority for now
When plotting the spectrum for the gravito_mhd
equilibrium the eigenvalues shift into the complex plane near the slow continuum point. This should not happen, all eigenvalues for this case must lie on the axes.
A similar issue occurs for the equilibrium defined by interchange_modes
.
One should naively expect that the eigenfunctions for the gravito-acoustic waves (HD) are all purely real, so that their imaginary parts are zero.
As of right now this is not the case (for some points the T1 eigenfunction has an imaginary part), so check what is happening here.
When both eigenfunction plotting and continua plotting are enabled, this sometimes raises a KeyError
when spectrum points are selected:
Traceback (most recent call last):
File "matplotlib/cbook/__init__.py", line 216, in process
func(*args, **kwargs)
File "legolas/python/visualisations/single_spectrum.py", line 59, in on_legend_pick
item = regions[artist]
KeyError: <matplotlib.lines.Line2D object at 0x118e033d0>
I suspect this has something to do with the event handler trying to activate/deactivate a figure legend based on the key corresponding to the spectrum point, throwing an error...
An error is thrown when not all eigenvalues have converged with shift-invert. It sometimes makes part of the spectrum unclickable or introduces a np.nan + 1j*0.0 solution at the location of the shift (like in this example).
Minimal example for reproduction
&gridlist
gridpoints = 100
/
&equilibriumlist
equilibrium_type = "RTI_theta_pinch"
/
&savelist
basename_datfile = "EP22_fig_3a30"
/
&solvelist
solver = "arnoldi"
arpack_mode = "shift-invert"
number_of_eigenvalues = 10
which_eigenvalues = "LM"
sigma = (0.49d0, 0.01d0)
/
Actual result
INFO | solving eigenvalue problem...
WARNING | ARPACK failed to converge! (maxiter reached)
WARNING | number of iterations: 100
WARNING | number of converged eigenvalues: 2 / 10
...
/Users/nicolas/codes/legolas/post_processing/pylbo/utilities/datfile_utils.py:478: RuntimeWarning: invalid value encountered in greater
eigenvalues[np.where(np.absolute(eigenvalues) > 1e15)] = np.nan
...
/Users/nicolas/codes/legolas/post_processing/pylbo/visualisation/spectra/spectrum_single.py:52: RuntimeWarning: invalid value encountered in greater
(self._nonzero_w_idxs,) = np.where(np.abs(dataset.eigenvalues) > 1e-12)
Version info
When running with the Arnoldi solver in some cases, I get an error:
ERROR | znaupd: starting vector is zero, try rerunning?
ERROR STOP
sometimes, re-running it indeed works as suggested, but usually it does not and the problem cannot be computed at all.
With no changes to the code, with a parameter file that reads:
&gridlist
gridpoints = 250
x_start = 0.0
x_end = 0.1
/
&equilibriumlist
equilibrium_type = "rayleigh_taylor"
/
&solvelist
solver = "arnoldi"
arpack_mode = "general"
number_of_eigenvalues = 10
which_eigenvalues = "LI"
maxiter = 16000
ncv = 400
/
&savelist
write_eigenfunctions = .true.
show_results = .true.
/
¶mlist
k2 = 1
k3 = 0
/
Actual result
INFO | the state vector is set to [rho, v1, v2, v3, T, a1, a2, a3]
WARNING | overriding x_end with 0.10000000
WARNING | standard equilibrium conditions not satisfied!
location of largest discrepancy (1): x = 0.02396385
value of largest discrepancy (1): 0.13546941E-11
amount of nodes not satisfying criterion (1): 984
INFO | ---------------------------------------------
<< Grid settings >>
geometry : Cartesian
grid start : 0.00000000
grid end : 0.10000000
gridpoints (base) : 250
gridpoints (Gauss) : 996
gridpoints (matrix): 4000
<< Equilibrium settings >>
equilibrium : rayleigh_taylor
wave number k2 : 1.00000000
wave number k3 : 0.00000000
default params : True
<< Physics settings >>
flow : True
external gravity : True
radiative cooling : False
thermal conduction : False
resistivity : False
viscosity : False
hall mhd : False
<< Solver settings >>
solver : arnoldi
ARPACK mode : general
number of eigenvalues : 10
which eigenvalues : LI
<< DataIO settings >>
datfile name : datfile
output folder : output
write matrices : False
write eigenfunctions : True
write derived eigenfunctions : False
write eigenfunction subset : False
---------------------------------------------
INFO | solving eigenvalue problem...
WARNING | Arnoldi: maxiter (16000) below recommended 10*N (40000)
ERROR | znaupd: starting vector is zero, try rerunning?
ERROR STOP
Matplotlib version 3.4.2 seems to lock zooming and panning on Pylbo plots after calling showall()
.
When doing this:
ds = pylbo.load("my_dataset")
p = pylbo.plot_spectrum(ds)
p.showall()
the pan and zoom buttons in the toolbar can be clicked, but zooming/panning does not work.
I'm assuming showall()
somehow disconnects some callbacks for these tools after the update.
Calling p.show()
instead of p.showall()
provides a workaround. To show all figures at once import matplotlib and call plt.show()
instead. The pylbo_wrapper.py
file calls showall()
, so this one needs to be modified.
In macOS Big Sur 11.6 it seems that even when we explicitly set write_eigenfunctions = .false.
the zgeev
routine in Lapack still references the right eigenvectors and hence throws a segmentation fault. This appears to be macOS-specific, Linux still works fine...
Maybe this is a bug in updated lapack/blas versions since everything worked fine before, so I'm gonna fix this with a dummy variable workaround.
Memory management seems extremely bad when datasets with eigenfunctions are loaded, without even loading the eigenfunctions themselves. This is probably a bug somewhere and should be tracked down. It makes no sense that 100s of MBs are used when no eigenfunctions are loaded into RAM.
Probably related to #45.
Example: running the discrete Alfvén modes with 201 gridpoints and eigenfunctions enabled:
import time
import pylbo
import memory_profiler
datfile = (pylbo.LEGOLAS_OUT / 'datfile.dat').resolve()
print('Memory before loading: {} Mb'.format(memory_profiler.memory_usage()))
t0 = time.time()
ds = pylbo.load(datfile)
print('Memory after dataset loading: {} Mb'.format(memory_profiler.memory_usage()))
t1 = time.time()
print('Time to load: {} sec'.format(t1 - t0))
efh = pylbo.EigenfunctionHandler(ds)
print('Memory after eigenfunction loading: {} Mb'.format(memory_profiler.memory_usage()))
t2 = time.time()
print('Time to load: {} sec'.format(t2-t1))
[output]
Memory before loading: [76.87890625] Mb
Memory after dataset loading: [746.9921875] Mb
Time to load: 0.8852338790893555 sec
Memory after eigenfunction loading: [1000.5546875] Mb
Time to load: 3.3105380535125732 sec
Memory after loading in the dataset jumps with 700 Mb, which should not happen.
I've noticed that a forced termination of a Legolas multirun with Pylbo does not always release the resources correctly. When forcing a KeyboardInterrupt
with ctrl-c
during program execution the subprocess correctly terminates, but the code seems stuck at trying to join the multiprocessing pool. This seems to happen independent of the amount of CPU's used.
This issue can be reproduced by calling pylbo.run_legolas()
for any configuration dictionary and immediately pressing ctrl-c
. In some cases termination is successful, in some cases it is not, so multiple attempts may be needed.
As reported by @jorishermans96, during generation of multiple parfiles with pylbo the creation of the basename_datfile
namelist item is not done correctly in some cases. If this key is not specified in the config dictionary supplied to pylbo.generate_parfiles
then pylbo keeps appending the full name.
Code for reproduction
import pylbo
config = {
"number_of_runs": 10,
"gridpoints": 50,
"parameters": {
"k2": 0.0,
"k3": np.pi * np.sin(0.25 * np.pi),
"cte_b03": np.linspace(0.001, 3, 10),
"cte_T0": 1.0,
"cte_B02": 0.0,
"cte_rho0": 0.0,
},
"equilibrium_type": "adiabatic_homo",
"show_results": False,
"write_eigenfunctions": True,
"write_matrices": False,
}
parfiles = pylbo.generate_parfiles(parfile_dict=config, basename="parfile")
The parfiles themselves are named correctly, but the basename_datfile
item in 0006parfile.par
looks like this:
basename_datfile = '000600050004000300020001parfile'
All tables in the radiative cooling / solar atmosphere are given in cgs and we simply normalise them.
Hence this is only correct for cgs units, and not for SI.
For the moment, using SI is locked via an error message to prevent erroneous use.
Answer tests should be regenerated when all modifications to investigate issue #13 are done.
Affected tests:
Others will be added here if affected
Look into optimisation of makefile instead of manually maintaining compilation order as is done right now. Take a look at cmake.
Compiling the code (release & develop builds) throws the following ranlib warnings in macOS
.../ranlib: file: ../lib/liblego.a(mod_version.f08.o) has no symbols
.../ranlib: file: ../lib/liblego.a(mod_physical_constants.f08.o) has no symbols
.../ranlib: file: ../lib/liblego.a(mod_version.f08.o) has no symbols
.../ranlib: file: ../lib/liblego.a(mod_physical_constants.f08.o) has no symbols
due to these modules not having procedures. Should be fixable by tweaking some CMAKE_RANLIB
settings in the build script somewhere.
When plotting quasimodes through the equilibria ideal_quasimodes
or resonant_absorption
the spectra look fine in a general sense, but the quasimodes do not show up.
This could be a resolution issue or equilibrium-related, not sure.
The default arnoldi solver mode is 'standard' on main and develop. However, on develop, this 'standard' mode has disappeared and hence the default should be changed.
Minimal example for reproduction
&gridlist
gridpoints = 10
/
&equilibriumlist
equilibrium_type = "adiabatic_homo"
/
&savelist
write_eigenfunctions = .true.
show_results = .true.
/
&physicslist
flow = .true.
/
&solvelist
solver = "arnoldi"
/
Actual result
INFO | solving eigenvalue problem...
ERROR | unknown mode for ARPACK: standard
ERROR STOP
Expected result
arpack_mode = "general" is picked by default
Version info
PR #69 added support for matplotlib 3.4, but apparently not all deprecationwarnings were properly handled. It seems that in some cases the warning is still triggered when the axes geometry is changed and this needs to be fixed before the 1.1 release.
Probably the easiest solution is to change how we add an additional subplot to the figure through _add_subplot_axes()
in the figure manager, this should propagate automatically to all instances that use it.
Running the adiabatic (or resistive) homogeneous case should give the same results in Cartesian and cylindrical geometries. However, in a cylindrical geometry the two outermost points in the spectrum are shifted further outwards.
Traced this back to mod_grid
, every value for x_start < 0.05
in cylindrical reproduces this issue. Quite positive the main issue occurs when setting the regularity conditions on r=0, since the subroutines defining the grids seem fine.
The only thing I can think of is that the matrix elements containing a factor 1/eps
blow up for small values, but this seems an unavoidable issue? The factor eps
will never be zero because it is given by grid_gauss
, so this should not matter much...
At some point we should set up code coverage in the CI framework.
Support is already there (building with cmake .. -DCoverage=ON
) so it will just be a matter of configuring the .yml file in Actions and linking to codecov or coverall.
The figures below shows the real part of the v1 and a2 eigenfunctions (not rescaled) for the first 5 eigenvalues in the slow mode sequence for the adiabatic homogeneous medium (Cartesian case).
The ones for a2 look fine, but for v1 there is something strange going on near the endpoints.
Pretty sure this is related (again) to the boundary conditions, and indirectly to #13.
One of the latest upstream changes to the arpack-ng
package changed the directory name
installed/include/arpack -> installed/include/arpack-ng
This apparently breaks our CMake link scripts, as cmake can no longer locate the arpack include directory.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.