Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

TwoPowerProfile #984

Merged
merged 12 commits into from
Apr 16, 2024
Merged

TwoPowerProfile #984

merged 12 commits into from
Apr 16, 2024

Conversation

ddudt
Copy link
Collaborator

@ddudt ddudt commented Apr 4, 2024

New class for profiles of the form $f(\rho) = a(1-\rho^b)^c$ where $a$, $b$, $c$ are all floats. This is equivalent to the "two_power" profile input format in VMEC.

@ddudt ddudt self-assigned this Apr 4, 2024
@ddudt ddudt added funtionality New feature or request to do things the code can't do now. EZ-review labels Apr 4, 2024
Copy link

codecov bot commented Apr 4, 2024

Codecov Report

Attention: Patch coverage is 94.59459% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 95.34%. Comparing base (d9acdfa) to head (bb6cae0).

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #984      +/-   ##
==========================================
- Coverage   95.34%   95.34%   -0.01%     
==========================================
  Files          87       87              
  Lines       21890    21926      +36     
==========================================
+ Hits        20871    20905      +34     
- Misses       1019     1021       +2     
Files Coverage Δ
desc/profiles.py 97.21% <94.59%> (-0.20%) ⬇️

desc/profiles.py Outdated
r = grid.nodes[:, 0]
if dr == 0:
f = a * (1 - r**b) ** c
# df/dr = inf at rho = 1 for all dr > 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

isnt this only true if b,c are floats <dr?

what is the use case for b,c being floats rather than ints?

my worry is also that if b<1 than the derivative blows up at r=0 as well.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True, depending on the parameter values it may or may not blow up at $\rho=0$ and/or $\rho=1$.

My personal use case is to use an integer b = 1 or 2, and c as a float. But if these numbers are optimizable params they all need to be floats.

Options:

  1. Add a check to make sure b >= 1. But there isn't much we can do about c.
  2. We don't need to merge this into the master branch. But I coded it up anyways so I figured I would share it.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My question is, do they all need to be optimizable? We could just make b,c constant integers and leave a optimizable?

Along a similar vein, VMEC has a bunch of other profiles types that might be worth looking at: https://github.com/jonathanschilling/educational_VMEC/blob/master/src/profile_functions.f

I have some VMEC input files that use the atan and two-lorentz profiles that would be good to try in DESC.

Copy link
Contributor

github-actions bot commented Apr 4, 2024

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     -2.47 +/- 3.51     | -4.60e-02 +/- 6.54e-02 |  1.82e+00 +/- 5.0e-02  |  1.86e+00 +/- 4.2e-02  |
 test_build_transform_fft_midres         |     -3.03 +/- 4.08     | -6.25e-02 +/- 8.43e-02 |  2.00e+00 +/- 5.3e-02  |  2.07e+00 +/- 6.6e-02  |
 test_build_transform_fft_highres        |     -1.94 +/- 3.25     | -4.77e-02 +/- 8.00e-02 |  2.41e+00 +/- 5.8e-02  |  2.46e+00 +/- 5.5e-02  |
 test_equilibrium_init_lowres            |     -2.03 +/- 3.93     | -2.05e-01 +/- 3.96e-01 |  9.87e+00 +/- 1.7e-01  |  1.01e+01 +/- 3.6e-01  |
 test_equilibrium_init_medres            |     -3.47 +/- 2.25     | -3.68e-01 +/- 2.38e-01 |  1.02e+01 +/- 1.6e-01  |  1.06e+01 +/- 1.8e-01  |
 test_equilibrium_init_highres           |     -1.68 +/- 2.16     | -2.09e-01 +/- 2.69e-01 |  1.22e+01 +/- 2.1e-01  |  1.24e+01 +/- 1.7e-01  |
 test_objective_compile_dshape_current   |     -2.13 +/- 10.49    | -7.84e-02 +/- 3.87e-01 |  3.61e+00 +/- 2.5e-01  |  3.69e+00 +/- 3.0e-01  |
 test_objective_compile_atf              |     -1.18 +/- 2.75     | -8.75e-02 +/- 2.04e-01 |  7.32e+00 +/- 1.6e-01  |  7.40e+00 +/- 1.3e-01  |
 test_objective_compute_dshape_current   |     -4.26 +/- 5.02     | -1.75e-04 +/- 2.06e-04 |  3.93e-03 +/- 7.7e-05  |  4.10e-03 +/- 1.9e-04  |
 test_objective_compute_atf              |     -0.20 +/- 7.27     | -3.65e-05 +/- 1.31e-03 |  1.80e-02 +/- 8.4e-04  |  1.80e-02 +/- 1.0e-03  |
 test_objective_jac_dshape_current       |     +0.24 +/- 9.42     | +1.03e-04 +/- 4.00e-03 |  4.26e-02 +/- 2.8e-03  |  4.25e-02 +/- 2.9e-03  |
 test_objective_jac_atf                  |     +0.27 +/- 3.52     | +5.09e-03 +/- 6.68e-02 |  1.91e+00 +/- 3.7e-02  |  1.90e+00 +/- 5.6e-02  |
 test_perturb_1                          |     +1.46 +/- 3.31     | +2.18e-01 +/- 4.96e-01 |  1.52e+01 +/- 2.3e-01  |  1.50e+01 +/- 4.4e-01  |
 test_perturb_2                          |     +1.52 +/- 4.56     | +3.04e-01 +/- 9.08e-01 |  2.02e+01 +/- 8.2e-01  |  1.99e+01 +/- 3.9e-01  |
 test_proximal_jac_atf                   |     +0.98 +/- 1.53     | +6.97e-02 +/- 1.10e-01 |  7.21e+00 +/- 9.6e-02  |  7.14e+00 +/- 5.2e-02  |
 test_proximal_freeb_compute             |     +1.98 +/- 0.82     | +2.52e-03 +/- 1.04e-03 |  1.30e-01 +/- 6.0e-04  |  1.28e-01 +/- 8.5e-04  |
 test_proximal_freeb_jac                 |     +0.32 +/- 1.61     | +2.33e-02 +/- 1.18e-01 |  7.37e+00 +/- 7.4e-02  |  7.35e+00 +/- 9.3e-02  |

desc/profiles.py Outdated
r = grid.nodes[:, 0]
if dr == 0:
f = a * (1 - r**b) ** c
# df/dr = inf at rho = 0 if b < 1
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is it even physical then if b < 1?

Copy link
Collaborator

@dpanici dpanici left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • add higher derivs
  • add to docstring that the params are only good for certain parameter ranges
  • add warning in build to warn if the parameters make things be divergent or nonphysical
  • check if it overlaps with a VMEC one and if so, then we can say it is similar to that

@ddudt ddudt changed the title PowerLawProfile TwoPowerProfile Apr 11, 2024
@ddudt
Copy link
Collaborator Author

ddudt commented Apr 11, 2024

  • add higher derivs
  • add to docstring that the params are only good for certain parameter ranges
  • add warning in build to warn if the parameters make things be divergent or nonphysical
  • check if it overlaps with a VMEC one and if so, then we can say it is similar to that
  • It turns out this is the same as the "two_power" profiles in VMEC, so I changed the name accordingly.
  • The MTanhProfile is also only implemented up to derivative order dr=2. Do we actually need higher order derivatives?
  • I added documentation and warnings for bad parameter values at initialization.

@ddudt ddudt requested review from f0uriest and dpanici April 11, 2024 20:46
@ddudt ddudt dismissed dpanici’s stale review April 12, 2024 13:59

made requested changes

f0uriest
f0uriest previously approved these changes Apr 13, 2024
Copy link
Collaborator

@dpanici dpanici left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see comments

desc/profiles.py Show resolved Hide resolved
desc/profiles.py Show resolved Hide resolved
@ddudt ddudt merged commit dc9d87f into master Apr 16, 2024
17 of 18 checks passed
@ddudt ddudt deleted the dd/power_law branch April 16, 2024 14:18
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
funtionality New feature or request to do things the code can't do now.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants