Code Monkey home page Code Monkey logo

legolas's Introduction

Hi there!

Twitter badge LinkedIn badge

Pinned repositories

Readme Card Readme Card

Stats

Niels's GitHub stats

Top Langs

legolas's People

Contributors

eprovst avatar jordidj avatar n-claes avatar nicolasbrughmans avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar

legolas's Issues

Bug in interactive plotting

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:

  1. Perform some filtering when selecting spectrum points, so that only points within a certain radius of the mouse are looked at. Possible issue: when the legend overlaps with the points.
  2. Make the legend a separate, standalone entity, such that clicking there is separated from clicking on the plot. This is probably the best course of action.

Low priority for now

Boundary conditions array derivative

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.

Documentation build is suddenly failing

Issue description

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...

Strange behaviour for interchange modes in multiruns

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.

Look into eigenfunction calculations/writing

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).

Minor discrepancy full resistivity/thermal conduction

There seems to be a small discrepancy in calculating the full resistivity and thermal conduction.

  1. Resistivity: it seems that we use a different prefactor pi instead of 4sqrt(2pi)/3, which differs by approximately 0.2. Also, log(coulomb_log) should actually be just coulomb_log.
  2. Thermal conduction: it looks like there is a discrepancy in converting between SI and cgs units for the perpendicular coefficient (the parallel one is fine).

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.

Pylbo: performance is bad for large datasets

At the moment, pylbo handles datasets as follows:

  • dataset information and basic data (eigenvalues, equilibria, parameters etc) are always loaded
  • eigenfunctions are (completely) loaded into RAM when queried
  • matrices are (completely) loaded into RAM when queried
    This works fine, however, larger datasets (> 600 mb) take a few seconds to load when trying to plot eigenfunctions and is quite an annoyance when iterating over results.

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:

  • make use of generators to handle eigenfunction loading
  • query eigenfunctions based on index, as we already do now
  • only load those queries into memory by iterating over the generator object

I have to test this in more detail, but I'm confident this will have a major impact on memory management and loading times.

Spectral webs

Extend the Python framework with the possibility of creating spectral webs of the results.
Before attempting this though it's probably a good idea to do these things first (in this order):

  • optimise the file storage format (#16)
  • interactive plotting bug (#20)
  • continua bug (#21)

setup script sometimes ignores user submodule

Issue description

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.

Bug report

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".

Additional cooling curve

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.

Coronal/photospheric flux tubes have wrong scaling

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

flux_tubes_legolas_roberts

Plotting of slow/Alfvén continua vs grid seems wrong

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...

Pylbo: inconsistent behaviour of `set_scaling`

Issue description

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.

Bug report

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.

spline functions

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.

Testing and documentation

Make sure that all submodules have extensive docs + regression tests, and that Pylbo has unit tests.
Current progress:

Docs

  • Legolas main routines documented
  • Legolas submodules documented
  • Pylbo code documented

Tests

  • Legolas main routines tested (core_tests)
  • Legolas submodules tested (regression)
  • Pylbo fully tested (pylbo_core)

Pylbo: cannot retrieve eigenfunctions in merged spectra if datfile contains only 1 eigenvalue

Issue description

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.

Bug report

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

  • Legolas version: 1.2.1
  • Pylbo version: 1.1.3
  • Python version: 3.8.13
  • Version of relevant libraries (matplotlib, numpy, etc.): matplotlib 3.5.3, numpy 1.23.3

Performance optimisations in solving linear systems

Issue description

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:

  1. zgbtrf to compute the LU factorisation
  2. zgbtrs to solve the linear system using the decomposition from 1.

As PR #97 introduced a lot of changes in the solver routines already, we should follow this order:

  • 1. Merge #97 into develop
  • 2. Rebase #98 onto develop and then merge
  • 3. Branch off develop and create the PR with this optimisation

That way the new changes can directly benefit from the updated regression framework.

Restructure equilibria

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.

Conduction & cooling equilibria

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.

Constant current multirun shows two additional branches

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...

Reduce unnecessary data usage

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:

  • Remove the left eigenvector matrix altogether and modify the QR and QZ routines accordingly
  • Modify the allocation routines to account for varying array sizes
  • Save size variables to the datfile, so Pylbo can use those for reference
  • Modify the Pylbo datfile reader to account for varying matrix sizes

All other stuff should not be affected, provided those subroutines use inferred arraysizes (which should be the case almost everywhere).

Segmentation fault for gcc <= 8

Issue description

For gcc <= 8 the current develop branch throws a segmentation fault at runtime (compilation is fine). From gcc 9 onwards everything runs fine.

Bug report

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.

Natural boundary conditions are wrong

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.

Optimise file storage format

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.

Eigenvalue selection from inverse vector-iteration spectrum

Issue description

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.

Bug report

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.
/

&paramlist
  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

  • Legolas version: develop branch, b7dbd98

Parallelised BLAS and LAPACK routines

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

Gravito-MHD has complex eigenvalues

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.

Gravito-acoustic (HD) eigenfunctions

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.

KeyError when plotting continua

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...

Shift-invert: error when not all eigenvalues have converged

Issue description

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).

Bug report

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

  • Legolas version: develop branch, b7dbd98

Error of "znaupd: starting vector is zero, try rerunning?"

Issue description

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.

Bug report

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.
/
&paramlist
  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 

pylbo: matplotlib 3.4.2 locks zoom after calling `showall()`

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.

Workaround

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.

segmentation fault when `write_eigenfunctions = .false.` in macOS

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.

Pylbo: memory issues

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.

Pylbo: resources are not released properly after forced termination

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.

pylbo: wrong datfile name generation for multiruns

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'

SI units are not consistent

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.

Restructure tests

Answer tests should be regenerated when all modifications to investigate issue #13 are done.
Affected tests:

  1. adiabatic homogeneous, due to 8ede7d3 (which changed equilibrium parameters as well)
  2. Suydam cluster modes, due to d978dba
  3. rotating plasma cylinder, due to d978dba
  4. kh_cd, due to d978dba
  5. rotating theta pinch, due to d978dba

Others will be added here if affected

Optimise makefile

Look into optimisation of makefile instead of manually maintaining compilation order as is done right now. Take a look at cmake.

Magnetothermal instabilities look quite wrong

First attempt at plotting magnetothermal instabilities, which don't look correct...
Eigenvalues are scaled as in right figure which comes from the paper (so the ones from legolas multiplied by -i) for a better comparison.

Screenshot 2020-03-19 at 14 23 03

cmake: no symbols warnings in macOS

Issue description

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.

Quasimodes do not show up

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.

Arnoldi solver default mode on develop

Issue description

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.

Bug report

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

  • Legolas version: develop branch, b7dbd98

Deprecationwarning: `change_geometry` in matplotlib 3.4+

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.

Cartesian-cylindrical at r=0

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...

CI code coverage

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.

Eigenfunction issue near the end points

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.

Screenshot 2020-03-12 at 16 25 59

Screenshot 2020-03-12 at 16 25 36

upstream ARPACK change broke CMake linking script

Issue description

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.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.