Comments (61)
To this comment, I want to add that I tried to run a fixed-boundary iter-like test case ( @lazersos gave me the input) that also does not converge, and by having a look at the interface geometries in the non-converged state, I also see that the innermost volume is going crazy, while the rest of the interfaces seem to be smooth. That lends support to Stuart hypothesis.
from spec.
@SRHudson wrote:
On geometrical regularization: see psifactor in preset.pdf and packxi.pdf. Start by reading the documentation, and as you have questions I will extend the documentation as required.
OK I will do that. Questions are coming soon.
@SRHudson wrote:
We need an input file (preferably not a free-boundary case, as there might be some other problems with free-boundary) for which spec fails to converge. This is usually because there is an interface close to the origin. An axisymmetric, stellarator-symmetric boundary might be ok, and would be faster for debugging etc.
Yes, I agree. I have a fixed-boundary case that does not converge (the "iter" case that Sam was once working on for rmp studies, which does not converge even in axisymmetry). Although this one is still non-stellarator-symmetric, shall I upload it in your favourite directory, namely TESTCASES?
In the meanwhile I can try to construct a tokamak toy-problem with an interface close to the axis, hoping we can make it fail to converge as well. Let me know.
from spec.
Can we change the name of the TESTCASES directory to something more ethical? How about TestCases, or even better InputFiles?
Joaquim: Please upload the non-stell.-sym. iter case; this will do for debugging. The lack of robust convergence is also associated with the angle constraints. Take a look at
http://w3.pppl.gov/~shudson/Spec/lforce.pdf
and then let's discuss. You should also learn how to compute the eigenvalues of the Hessian, see
http://w3.pppl.gov/~shudson/Spec/hesian.pdf
from spec.
Can we change the name of the TESTCASES directory to something more ethical? How about TestCases, or even better InputFiles?
Done.
Please upload the non-stell.-sym. iter case; this will do for debugging.
Done.
Take a look at http://w3.pppl.gov/~shudson/Spec/lforce.pdf and then let's discuss.
Two questions:
(1) I understand how the coordinate axis is calculated, and I understand that the extrapolation from the innermost interface towards the coordinate axis is done in such a way the the high-m modes decay rapidly towards the axis, while the m=1 mode goes like r~sqrt(s). But are the high-m modes decaying fast enough? Also, why are the m=0 (n>0) modes decaying like ~ s, and not another power?
(2) Having a look at the papers by Hirshman et al you refer to in the documentation of lforce, I can see that the spectral condensation is supposed to be implemented in your spectral constraints as the second term in Eqs. (3) and (17) of lforce.pdf. But what are the other terms (meaning and derivation) and why are they clustered together with the spectral constraint? Unless the coefficients alpha,beta,gamma,delta are Lagrangian coefficients, but then they would not be user-supplied but output values right?
from spec.
- I pulled the master and found the input files, thanks. Are the input files available on other branches?
- There are three regularizations near the coordinate axis. The first is how the interface geometry is extrapolated into the coordinate origin, see
http://w3.pppl.gov/~shudson/Spec/coords.pdf
the second is the inclusion of a "pre-conditioning factor", see
http://w3.pppl.gov/~shudson/Spec/packxi.pdf
and the third is the inclusion of regularization factors on the vector potential, see
http://w3.pppl.gov/~shudson/Spec/matrix.pdf
Joaquim: let's revise this together. Let's start with the interpolation of the interface geometry into the coordinate axis. Let me ask you: is the coordinate label, s, used in spec in the innermost volume proportional to the "minor radius" of a coordinate surface, or is s proportional to the toroidal flux? You might like to plot surfaces of constant s for the simple model r(s,t) = s cos(t), z(s,t) = s sin(t), and for r(s,t) = \sqrt s cos(t), z(s,t) = \sqrt sin(t). (We can discuss how the spectral condensation constraint does not constrain the location of \t = 0 later.)
from spec.
- I corrected what I think was an error in the iterfixed.sp file, and now it converges. Note that since
"Existence of three-dimensional ideal-MHD equilibria with current sheets"
J. Loizu, S.R. Hudson et al., Phys. Plasmas 22, 090704 (2015)
SPEC has allowed for discontinuous rotational-transform profiles, and this is accomplished by including two values for the rotational-transform, iota and oita, for each interface, one for each "side" so to speak. In iterfixed.sp, oita was not provided; so the input rotational-transform profile was weird. So the first task was to set oita = iota. Second, I had to play around a little, with Linitialize=1, Lfindzero=0 and Lconstraint=1, to find suitable values for mu and pflux, the helicity multiplier and the enclosed poloidal flux in each volume. This can be a little tricky, as there are a discrete family of solutions . . . Then, we can set Lfindzero=2 to iterate on the geometry to find global force-balance. I have pushed the new copy of iterfixed.sp to the master branch. Please pull and check everything.
from spec.
I understand the issues with the input, thanks! I pulled, compiled and run the master branch on the fixed-boundary iter case and indeed convergence is successful. Attached a poincare plot with the interfaces and the iter vessel for illustration. Sam can now use this as a basis for rmp studies.
from spec.
Let's start with the interpolation of the interface geometry into the coordinate axis. Let me ask you: is the coordinate label, s, used in spec in the innermost volume proportional to the "minor radius" of a coordinate surface, or is s proportional to the toroidal flux? You might like to plot surfaces of constant s for the simple model r(s,t) = s cos(t), z(s,t) = s sin(t), and for r(s,t) = \sqrt s cos(t), z(s,t) = \sqrt sin(t). (We can discuss how the spectral condensation constraint does not constrain the location of \t = 0 later.)
According to Eqs.(7) and (8) of
http://w3.pppl.gov/~shudson/Spec/coords.pdf
the "minor radius", r, of a coordinate surface is, as I said before, r~sqrt(s) since that is how the m=1 modes decay towards the axis. Therefore, the label s is proportional to the toroidal flux, s ~ r^2, at least to lowest order (i.e., neglecting m>1 modes). What I do not get so much is why the m=0 modes have a regularization factor that goes linear with s. Otherwise, I think this regularization makes sense.
Regarding the "pre-conditioning factor",
http://w3.pppl.gov/~shudson/Spec/packxi.pdf ,
why is it necessary and how exactly is it determined? I see that, in preset.h,
preset.h:935: psifactor(1:mn,1:Nvol) = one
preset.h:941: if( im(ii).eq.0 ) then ; psifactor(ii,vvol) = tflux(vvol)( +half) ! 28 Jan 15;
preset.h:942: else ; psifactor(ii,vvol) = tflux(vvol)(halfmm(ii)-half) ! 28 Jan 15;
preset.h:952: if( im(ii).eq.0 ) then ; psifactor(ii,vvol) = Rscale * tflux(vvol)**zero ! 08 Feb 16;
preset.h:953: else ; psifactor(ii,vvol) = Rscale * tflux(vvol)**halfmm(ii) ! 29 Apr 15;
from spec.
Try interpolating -- any way you like -- between the coordinate axis
R = 0.95, Z=0.00,
and the interface
R=1.05 + 0.3 cos t + 0.1 cos 2t + 0.05 cos 3t, Z = 0.3 sin t
or something like that. I just guessed these numbers, you should play around with some different interfaces, different coordinate axes including non-stellarator-symmetric. The coordinate interpolation is somewhat arbitrary: all that is required is that the coordinate surfaces do not intersect.
I just guessed the pre-conditioning factor. I cannot provide any conclusive argument why it is required and what form it should take. The only way to determine this is to examine whether the global convergence improves or not with this regularization factor. Note that the regularization factor and the coordinate interpolation into the coordinate axis are completely separate.
from spec.
I agree that s ~ r^2. This has implications for the regularization factors on the vector potential.
Regarding psifactor, we can easily include an input option, e.g. ifactor = 0, 1, 2, . . . that would allow the user to select from a variety of choices. Then, we could systematically explore any benefits or disadvantages. When I included this factor, spec seemed to converge more robustly, but I did not perform an exhaustive investigation.
Perhaps you merge ma00aa into master, and then delete this branch. Then, create another branch, e.g. psifactor, and I will include an input flag etc. so that you can investigate this more thoroughly.
from spec.
Perhaps you merge ma00aa into master, and then delete this branch. Then, create another branch, e.g. psifactor, and I will include an input flag etc. so that you can investigate this more thoroughly.
Before I merge the ma00aa into master, I would like to see some feedback on the issue "comparing two output files", since the two branches are not producing the exact same result for a given input.
from spec.
Try interpolating -- any way you like -- between the coordinate axis
R = 0.95, Z=0.00,
and the interface
R=1.05 + 0.3 cos t + 0.1 cos 2t + 0.05 cos 3t, Z = 0.3 sin t
or something like that. I just guessed these numbers, you should play around with some different interfaces, different coordinate axes including non-stellarator-symmetric. The coordinate interpolation is somewhat arbitrary: all that is required is that the coordinate surfaces do not intersect.
Well, obviously if the interface geometry is really crazy and the interpolation does not damp high-m modes rapidly enough, then surfaces can overlap. Playing with examples is useful to get a sense, but will not lead to finding the best or robust interpolation method.
There is, however, a general "sine qua non" condition, and that is that coordinate surfaces will not overlap if and only if the Jacobian of the coodinates is positive everywhere.
Given that the Jacobian is given by Eq.(18) of
http://w3.pppl.gov/~shudson/Spec/coords.pdf
this condition translates into a general condition that
(dZ/ds)(dR/dtheta) > (dR/ds)(dZ/dtheta)
but I do not see any way of deriving an explicit condition for the interpolation power law, even if one plugs in Eqs.(7) and (8) in this expression.
Nevertheless, I believe that the interpolation method you implemented shall work in most cases, and if that would be the problem, then wouldn't that appear to us clearly? Since sqrt(g)<0, namely negative jacobian, would produce an error?
from spec.
Regarding psifactor, we can easily include an input option, e.g. ifactor = 0, 1, 2, . . . that would allow the user to select from a variety of choices. Then, we could systematically explore any benefits or disadvantages. When I included this factor, spec seemed to converge more robustly, but I did not perform an exhaustive investigation.
We can try to explore this, but I would like to understand why this came to your mind at all? For example, does VMEC also apply such kind of factor?
from spec.
Interpolation
If you can derive an interpolation procedure that guarantees that the interpolated surfaces do not overlap then that would be fantastic.
I know a method that will work, but that will be impractical. Consider the solution to Laplace's equation, \div\grad\phi=0, with the boundary condition \phi=0 on the coordinate axis (or a toroidal surface very close to the coordinate axis) and \phi=1 on the innermost interface. The surfaces of constant \phi are guaranteed to not-intersect. Does anyone have a coordinate-independent subroutine for solving Laplace's equation in 3D? I tested this idea in 2D using
http://www.nag.co.uk/numeric/FL/manual19/pdf/D03/d03eaf_fl19.pdf
and it works beautifully. However, as I mentioned, I think that this approach will be too slow in SPEC. Ideally, SPEC could compute the Beltrami fields using integral equation method, which don't require coordinates! Antoine Cerfon NYU has published various papers on this topic. I don't think he has gone "three-dimensional" yet. (In fact, it might be a good idea to invite Antoine to join this Git project. What do you think?)
pre-conditioning
Imagine making a change of size 1.0e-02 to an m=0 or an m=1 mode. Not such a crazy thing to do, right? Imagine making a change of size 1.0e-02 to an m=24 mode. Something about this seems crazy, as the surface now becomes extremely rippled, and this would cause lots of problems. The "regularization" psifactor was intended to avoid this.
I don't know what VMEC does, but I would be interested to learn.
The only criteria of concern is, does this factor make SPEC more robust or not? We can play around with this factor and see what works, include a choice of options for the user, or theoretically derive the best method. We want the hessian matrix to be as close as possible to the identity matrix. So, take a simple circular cross section tokamak with lots of surfaces, compute the equilibrium, and then examine the eigenvalues of the hessian and see what the effect of this factor is. This will be a very fruitful exercise that is almost guaranteed to improve the convergence of SPEC. This analysis will also be a major step forward in extending SPEC to calculate linear stability. Note that the hessian matrix is already being computed by SPEC, and note that this matrix presently also includes the spectral constraints. In fact, the only way to understand the convergence properties of SPEC is by understanding the hessian. Also, consider this: imagine that we have computed an iter equilibrium, for example, and we want to make a small change (in the pressure, rmp, etc.). Well, a fantastic pre-conditioning matrix will be the inverse hessian!
from spec.
Interpolation
I understand that a coordinate-independent solution of Laplace's equation in 3D would solve the problem. But since this coordinate grid is generated at each iteration on the geometry, it seems that this would make things very slow, as you pointed out. But of course bringing Antoine Cerfon could be interesting at some point.
I would say for the time being let us identify which part in the code makes the convergence fail or non-robust. I am going to have a look at the non-converged case of free-boundary iter and see if the innermost coordinate geometry looks good or not. If it still looks good despite the non-convergence, I would reckon that the problem lies somewhere else, e.g. the preconditioning.
pre-conditioning
I am now reading the paper by Hirshman on the preconditioning in vmec,
Hirshman_1991.pdf
which I guess is based on the same idea?
I am also having a look at the spec source in order to understand better what exactly does the code do with this psifactor.
from spec.
Interpolation
Have a look at the grid in the non-converged iterfree case (i.e. when Linitialize=1), here shown for the innermost volumes,
iterfree_poincare_grid.pdf
In pink are the grid points, in red the interfaces, in black poincare data points. It looks like eventhough the innermost interface goes crazy, the coordinate surfaces remain quite smooth, i.e. no crossing.
from spec.
The idea of a coordinate-independent Laplace solver was not so much for the construction of the coordinate grid, but if we can construct a coordinate-free method for Laplace in 3D then we could probably have a coordinate-free method for constructing the Beltrami fields directly, and not bother about the coordinates at all. Furthermore, we usually only need B^2 on the interfaces, so it is not required to calculate the magnetic field globally (as we are doing now).
from spec.
The idea of a coordinate-independent Laplace solver was not so much for the construction of the coordinate grid, but if we can construct a coordinate-free method for Laplace in 3D then we could probably have a coordinate-free method for constructing the Beltrami fields directly, and not bother about the coordinates at all.
I see, but the beltrami field (more precisely the corresponding vector potential) only satisfies a Laplace equation if mu=0 (i.e., in vacuum) right?
from spec.
Furthermore, we usually only need B^2 on the interfaces, so it is not required to calculate the magnetic field globally (as we are doing now).
But usually to get B^2 on the inner side of the innermost volume, one solves for B in the volume first. Is it possible to get B^2 on that side without solving for B inside the volume?
from spec.
The Beltrami equation and Laplaces equation are not too much different, so I thought that if we understood now to implement a coordinate-free approach to one then we should be able to understand how to implement a coordinate-free approach to the other. This is what Cerfon is working on. With such coordinate-free, Green's function approaches, it is not required to obtain the solution globally, so as I understand we could solve for B only where it is required. The fact that a coordinate grid is not required and we need only compute B where we need makes this approach very attractive. This is why Cerfon's work could be a great benefit to SPEC. Note that if we need to compute the helicity integral, then we will need both B and A globally.
from spec.
I would like to make my point very clear: the problems with convergence, whether the spectral constraints are appropriate and implemented correctly, and understanding of the "pre-conditioning" factor can all be understood by examining the eigenvalues and eigenvectors of the hessian. The hessian is already being computed by SPEC.
Imagine turning all the spectral constraints off. Then the tangential degrees of freedom in the Fourier representation should correspond to eigenvectors of the hessian with zero eigenvalues. Then we can add the spectral constraints back in one by on, and see the zero eigenvalues disappear. Note the the spectral condensation algorithm of Hirshman et al. as is implemented in SPEC does not uniquely constrain the Fourier representation (the location of the \t = 0 line is not constrained) and this leads to convergence issues. There are a couple of other issues aswell (what if the boundary is not given in the spectrally condensed angle?), which is why I included other constraints, see
http://w3.pppl.gov/~shudson/Spec/lforce.pdf
All of the coding (for the hessian, spectral constraints) has been implemented, but a careful, rigorous revision and perhaps debugging is required. The answers are all to be found in the eigenvalues and eigenvectors of the hessian! This information will also provide information about the linear stability of SPEC equilibria, and extracting that information would be a game changer.
from spec.
the problems with convergence, whether the spectral constraints are appropriate and implemented correctly, and understanding of the "pre-conditioning" factor can all be understood by examining the eigenvalues and eigenvectors of the hessian.
Well then, let's calculate the eigenvalues and eigenvectors of the hessian! Shall I work on this?
Also, shall we construct a simple toy-problem (circular tokamak cross-section, for example) to debug and analyze the hessian? I just worry that if the geometry is too simple, there will be no convergence problems.
Note the the spectral condensation algorithm of Hirshman et al. as is implemented in SPEC does not uniquely constrain the Fourier representation (the location of the \t = 0 line is not constrained) and this leads to convergence issues.
Does this mean that these issues should also be present in VMEC?
from spec.
Yes, Joaquim, because you have time to concentrate on this, I would like you to examine the hessian, and how the various spectral constraints and pre-conditiong factor affect the eigenvalues/vectors. I suggest (1) construct a simple circular cross-section equilibrium (including appropriate spectral constraints); (2) keeping the equilibrium fixed, calculate the hessian with and without both the preconditioning factor (i.e., set psifactor=1) and the spectral constraints; (3) and investigate . . . We don't want need to worry about whether the code converges or not because we are examining the hessian directly.
I don't know enough about VMEC to make any conclusive comments about how the spectral constraints are implemented, but I did at one stage convince myself that, because SPEC has a slightly different implementation to that of VMEC, that there were some issues that SPEC needed to address that was not present in VMEC.
from spec.
(1) construct a simple circular cross-section equilibrium (including appropriate spectral constraints)
Done. I uploaded an input file to the InputFiles/TestCases/ directory (on the master branch). Below I attach a poincare plot of the converged equilibrium.
from spec.
(2) keeping the equilibrium fixed, calculate the hessian with and without both the preconditioning factor (i.e., set psifactor=1) and the spectral constraints;
The Hessian is computed by SPEC and written in the .DF file, if I understand correctly. Also, there are two flags, according to your documentation, namely LHevalues and LHevectors, that can be set to 'true' in order to calculate the eigenvalues and eigenvector of the Hessian. I tried this, and got the error message:
global : fatal : myid= 3 ; .true. ; eigenvalue solver needs updating to F08NAF ;
which originates because there is, in hesian.f, a check
#ifdef NAG18
call F02EBF(...)
#else
error message
So, shall I try to update the source in order to call F08NAF?
from spec.
Joaquim:
- Construct a simple equilibrium in Cartesian geometry and study the hessian. There are no spectral contraints, and no geometry so-to-speak. Introduce a preconditioning factor that reduces the hessian to the identity matrix.
- Construct a simple equilibrium cylindrical geometry, I have not run spec in in cylindrical geometry for a while, but I guess that all the subroutine will still work. Examine the hessian, and introduce a preconditioning factor that reduces the hessian to the identity.
- Construct a simple equilibrium in toroidal geometry; the simplest non-trivial case that you can (i.e., circular cross-section). Study the hessian, and introduce a preconditioning factor that reduces the hessian to the identify. (In the toroidal case, the preconditioning factor will be a matrix. In the Cartesian and cylindrical case, the preconditioning matrix will be diagonal.) In the large aspect ratio limit, this preconditioning matrix should reduce to the cylindrical preconditioning matrix. (We will need to understand everything about the spectral constraints, which only appear in toroidal geometry.)
- Usually in toroidal geometry, the preconditioning matrix for a circular cross-section equilibrium will be sufficient for arbitrary toroidal geometry, and this could be used as the default for spec. I think that we could determine an analytic approximation to this preconditioning matrix. I would just fit some analytic form to the numerically constructed inverse hessian, or we could try to derive something from first principles like Betancourt did.
- The most powerful will be to implement the following. Imagine that an axisymmetric e.g. iter-like equilbrium has been constructed, or similarly a W7-X equilibrium. The best preconditioning matrix will be the inverse of the hessian, which spec computes. So, save the inverse hessian to file, and for all calculations with a similar geometry we can read this "best" preconditioning matrix from file.
from spec.
F08NAF?
Please either see what you have available in your NAG library, or mention to Josh that this NAG replacement is a priority. Eventually we intend to eliminate NAG, so to reduce effort it is probably best to find an alternative to NAG immediately. The hessian is symmetric I believe (please check) and this can be exploited. Futhermore, the hessian is block tridiagonal. Exploiting this will lead to a significant reduction in cost.
I also think that we should merge the NAG replacement branch with the master branch sooner rather than later. I can't help but think that this will reduce future conflicts, and it helps the NAG replacement activity if we have the faster spec (from the ma00aa branch) and that we have all the "fixes". For example, Josh had to re-discover that the dspec needs the -r8 flag in the Makefile. This was fixed in the master, but not the NAG_REPLACE branch. We should not have to fix the same bug twice.
As for your tokamak equilibrium, I think that we could identify the generic structure of hessian if the equilibrium were as close as possible to having circular flux surfaces (i.e., perhaps we should start with an equilibrium without pressure and at large aspect ratio).
from spec.
F08NAF?
yes, that is what you wrote as a message: "update to F08NAF"
Please either see what you have available in your NAG library
I have F08NAF in my library
Eventually we intend to eliminate NAG, so to reduce effort it is probably best to find an alternative to NAG immediately.
True. But I have no idea of the time-scale in which NAG will be replaced - especially since it seems that time availability is a concern for most of you. So, perhaps I can, for the time being, just use F08NAF. I just need to learn how to call it.
I also think that we should merge the NAG replacement branch with the master branch sooner rather than later.
Let us assume we merge (after testing) now the nag_replace branch, then what? Any further "nag-replacement-improvements" should be done on a separate branch than master, so we would not delete the nag_replace branch. Thus the only thing we would win is that master would have some nag routines replaced. Not sure this is useful. Let me know what do you think.
As for your tokamak equilibrium, I think that we could identify the generic structure of hessian if the equilibrium were as close as possible to having circular flux surfaces (i.e., perhaps we should start with an equilibrium without pressure and at large aspect ratio).
I see what you mean. Although apparently you are suggesting starting from a slab case, so I better do that rightaway.
from spec.
Just a comment. @jloizu If you are trying to find a eigenvalues solver. You can try F08FAF, which I use in FOCUS.
And for similar jobs, what I did is to write the Hessian into a file (I used hdf5) and then use matlab or python to do eigenvalue decomposition and analysis stuffs.
I'm really looking forward to seeing your results.
from spec.
Caoxiang raises a point worth considering. The eigevalues etc. of the hessian are, at present, purely used for diagnostics only, and so it would simplify the spec source if this were "out-sourced" to an external routine. However, looking ahead, we plan to develop spec into a linear stability code, and in this case it would be good to keep the eigen-analysis internal to spec. Spec would then compute both the equilibrium and the linear stability.
from spec.
As a compromise I decided to write my own matlab routine to the read the hessian binary file (.DF) and extract the matrix, then I can analyze it easily with matlab. In the near future, when including such calculations in SPEC, we may directly avoid NAG.
from spec.
For the circular tokamak equilibrium, here is what I get for the Hessian matrix.
First attachment shows the amplitude (in absolute value and in log-scale) of the matrix elements. Clearly the matrix is essentially block-tridiagonal. However, it seems that the values of the matrix elements increase as we go "down" the diagonal, which is the opposite I expected (although this may be related to the presence of the spectral contraints). Notice that the "matrix" is shown such that the (1,1) coefficient is on the bottom left corner.
Second attachment shows the actual values of the matrix elements, but I saturate the colorbar at 0.01 in order to see structure. Actual maximum values go up to +/- 1e4. Strangely, there is an oscillatory structure all along the tridiagonal (values jumping from positive to negative very rapidly).
I guess analyzing a slab case will help geting a better understanding.
from spec.
Another remarkable result is the same Hessian plots obtained for the same equilibrium but with pscale=0, i.e. with no pressure (the shafranov shift is then smaller and surfaces more circular).
circulartok_nopressure_hessian1.pdf
circulartok_nopressure_hessian2.pdf
The absolute values and structure of the diagonal remains similar, but the Hessian is no longer tridiagonal!
from spec.
That's interesting. I looked at my results of FOCUS. The Hessian matrix is beautifully dominant block-tridiagonal.
This is an result of 16 circular coils producing circular torus, without spectral condensation. The 16 coils are all the same during the computation. That's why the block size is 16.
And I also have a result of same inputs but turning on the spectral condensation.
Roughly, they look similar. I should go further. But I'm distracted to finish something else.
It's the result from my code FOCUS. Just for your reference.
from spec.
- I suggest you simplify the equilibrium as much as possible. Perhaps you can use only one internal interface and low Fourier resolution.
- After constructing the equilibrium, calculate the hessian with all the spectral constraints turned off. This might be possible by setting opsilon, epsilon and upsilon each equal to zero. You should become familiar with looking at the source. I think that lforce.h is easy enough to read (though you will need to read through carefully), but dforce.h is messy. In dforce.h, look for dFFdRZ and hessian. I hope that you can see what is happening.
- Before we analyze the hessian too much, we should make sure that it is correct. One test is to compute the hessian using finite-differences. See Lcheck=5. Also, to make things simple, we can compute the hessian under the constraints of fixed poloidal flux and fixed helicity multiplier, i.e. for Lconstraint=0. With Lconstraint=1 it is required to compute how these two parameters must vary as the geometry is varied in order to satisfy the rotational-transform constraint, and this is an extra complication.
- For each eigenvalue, please construct a file showing the eigenvector, e.g. {\bf v} = (R_{m,n},Z_{m,n})^T. With the spectral constraints turned off, there should be a tangential null space (but be careful, the null space might only be resolved at high resolution). Try to identify each eigenvector in terms of a deformation in real space.
- Start with Cartesian. There are no spectral constraints, the geometric deformations are easy to identify, . . .
- Just investigate everything that you can!
from spec.
- I found a bug in preset that I introduced recently. I used lp as a loop variable, but lp is also an input variable. This bug might only affect cases for which the rotational transform, oita, is specified by the lp, lq, rp and rq arrays. I have corrected the bug and pushed the fix to master.
- I have added an input file in Cartesian geometry.
- I un-commented out the eigenvalue routine, F02EBF, which seems to work for me, in hesian.h. If this doesn't work for you, just comment it out and perform your own eigenvalue analysis using matlab.
- For the input Cartesian file, I have turned on the input flags to check the analytic calculation of the hessian, and it looks ok (in that the finite difference derivatives seem to agree with the analytic values).
- I computed the eigenvalues and eigenvectors for this Cartesian case. The eigenvectors are simple: v_m = R_m. There is no toroidal or poloidal coupling in this very simple case. The eigenvalues seem to satisfy \sqrt(\lambda_m) = a + b m for m>0. Please investigate this for yourself.
- I changed psifactor in preset (search for "08 Aug 17"). The eigenvalues now satisfy \lambda_m = const. for m>1. Isn't this what preconditioning is all about?
- I have committed all these changes to master. Sorry, but in this case there was a good reason. First, there was a bug, and I am very confident that I fixed this bug. Second, I commented out the changes I made to preset, so there was in fact no change. Third, the changes I made to hesian were trivial screen output. I will in future create a new branch whenever I make a change so that you can review.
- We need lots more example input files for testing . . . and we need a script that will automatically run spec on all these input files. In the above I described how I recently introduced a bug that did not have any impact on the l2stell or the tokamak cases. Ideally we would have enough input files that would exercise and check every capability, every subroutine of spec. We should also have a standardized way of naming the input files.
- Please check that you can reproduce all of the above statements regarding the eigenvalues etc. and how changing psifactor affects the eigenvalue spectrum. Then, explore the effect of including more interfaces, etc. etc. I am just guessing my way through, and I am leaving it to you to check that everything is done precisely and rigorously.
- If we can find an unstable Cartesian equilibrium then we can write a great paper. Imagine constructing a sequence of axisymmetric equilibria with increasing pressure, for example. For each equilibria, we can re-run spec allowing for non-axisymmetric harmonics, and if we can find an axisymmetric equilibrium that gives negative eigenvalues of the hessian, then we have demonstrated spec's capability as a linear stability code for partially relaxed equilibria (i.e., equilibria with islands).
from spec.
I un-commented out the eigenvalue routine, F02EBF, which seems to work for me, in hesian.h. If this doesn't work for you, just comment it out and perform your own eigenvalue analysis using matlab.
This is not how we should work. This is bad programing philosophy. I pulled the master branch, and now the source does not compile.
hesian_r.o: In function
hesian_': hesian.F90:(.text+0x33ef): undefined reference to
f02ebf_'
simply because the F02EBF nag routine is obsolete in NAG Mark 24, which is the version of NAG we use here. Of course this is an easy fix, e.g. un-comment the #ifdef that you commented.
But we don't want that each of us has a different version of SPEC. That is why the master branch should ALWAYS be compiling and running for all users, and that is why we should create short-lived branches that can be first tested and then merged. The only cases where applying changes directly to master make sense are when adding an input file, or something equivalent. This is my philosophy. If you have another one, then persuade me that it is better.
from spec.
I fixed this problem by doing
git checkout 075c256 hesian.h
git commit -m "Reverted the file hesian.h to..."
git push origin master
which reverts (only) the file hesian.h to the last commit before that change was made. Also, the history is not modified, so every change that was made is still traceable.
@SRHudson if you want to make a change in hesian.h such that F02EBF is called, then maybe define the environment library NAG18 and add it as compilation option (only for pppl compiling option), something like "mpif90 -DNAG18 ..." should work. Or whatever you want, but let's create a branch and test always before merging! :)
from spec.
I am working on the cartesian case. @SRHudson, I see that in the input file you uploaded,
forcetol=1e+12
which I guess is a typo (not important if the system is at equilibrium already though, I suppose)
from spec.
I did run the cartesian case you provided, setting forcetol=1e-12, and running once with Lcheck=0 to write the .DF file, then with Lcheck=5 to use your checking routine.
The Hessian looks good, the diagonal seems to agree also with your analytical estimate given by Lcheck=5, but there are two corners, namely H(1,9) and H(9,1), which differ substantially. Below what I get for the Hessian matrix elements from the SPEC binary output file .DF,
(now the plot really shows the matrix in the conventional form, i.e. H(1,1) is the upper-left corner)
By definition, this matrix should be symmetric, right?
from spec.
Looking at dforce.h, I see that the hessian is assigned as:
(line 1449)
hessian(tdoc+idoc+1:tdoc+idoc+LGdof,tdof) = - dFFdRZ(idoc+1:idoc+LGdof,vvol+0,1,idof,0)
where
(line 1037)
dFFdRZ(idoc+1:idoc+mn ,vvol,iocons,idof,innout) = + efmn(1:mn ) * psifactor(ii,vvol-1+innout) * BBweight(1:mn)
and I am trying to see why I do not get a symmetric matrix. Some parameters I do not understand, what is for example LGdof? How is it different from mn?
from spec.
"This is not how we should work. This is bad programing philosophy."
I am sorry. I agree with your comments. I was lazy and arrogant, and I will not let this happen again.
I will address your other comments/questions as soon as I can (perhaps days). forcetol = + huge means that spec will be satisfied that the provided input is in equilibrium and will make no effort to change the geometry. LGdof = local geometric degrees of freedom to each interface. Each interface has an R and Z, and for stellarator symmetric the R function has mn harmonics and the Z function has mn-1 harmonics (the sin(0\t-0\z) is irrelevant). Let me ask you: where in the code is LGdof assigned? What is NGdof? Where is mn assigned? Clue: these are global variables that only need to be assigned once (perhaps preset is a good place to look, but some global variables need to be set in global:readin before the entire input file is read from file).
I think that the hessian should be block tridiagonal, and for this Cartesian case (where there is no toroidal/poloidal coupling) each block should be diagonal. Please investigate. Packing and unpacking the geometrical degrees of freedom is required to allow standard numerical routines to be employed, and this "packing" should be understood. See
http://w3.pppl.gov/~shudson/Spec/packxi.pdf
from spec.
Thanks for the initial comments.
I think that the hessian should be block tridiagonal, and for this Cartesian case (where there is no toroidal/poloidal coupling) each block should be diagonal.
I would like to understand first, regardless of the numerics, what is the definition of the Hessian matrix?
Is it defined as
H_{ij} = (d^2 F)/(dx_i*dx_j)
where F is the energy functional and {x_i, i=1...NGdof} are the geometrical independent variables? If that is the case, then the matrix should symmetric.
Or is it defined as
H_{ij} = (d f_i)/(dx_j)
where f_j labels the harmonics of the force imbalance? If that is the case, I am not so sure that the matrix should be symmetric.
from spec.
I will investigate as you suggest. I just want to let you (all) know that I will be on holiday for 3 weeks, August 11th - August 31st. I will use all my time on this when I am back.
from spec.
I discovered something.
To illustrate it, consider taking the cartesian.sp test case uploaded by @SRHudson, and make this changes before running:
mpol = 22
forcetol = 1e-12
lhevalues = F
lhevalues = F
lcheck = 0
linitialize = 1
c05xtol = 1e-6
Running this case produces a converged equilibrium after a few iterations, with only one Hessian evaluation required. In that case, the hessian is perfectly diagonal, and with eigenvalues increasing with m, as @SRHudson observed previously.
cartesian_hessian_xtol1em6.pdf
Now, delete the Hessian file, "rm .cartesian.sp.DF", and rerun the same case but with
c05xtol = 1e-10
and repeat the same procedure for
c05xtol = 1e-12
which produces two other converged states (with lower force imbalance), after 2 and 3 Hessian evaluations, respectively. In this case, however, the Hessian is no longer diagonal, nor symmetric, even though the diagonal terms remain similar.
cartesian_hessian_xtol1em10.pdf
cartesian_hessian_xtol1em12.pdf
The way I understand this is as follows: in the first case the Hessian is evaluated w.r.t an equilibrium that is perfectly symmetric (only m=0,n=0 modes) and thus there is no coupling possible between modes, thus Hessian is diagonal. However, as force is brought down to zero, perhaps some modes become, even if very small, finite, and thus re-evaluation of the hessian on a state that is not exactly symmetric will produce coupling. Thus the Hessian starts becoming non-sparse.
Not entirely sure this is expected/physical, but it is there. And even if it is expected, that is important information because the pre-conditioning must take into account that the Hessian is not diagonal or block-tridiagonal.
Hope this gives some insights.
from spec.
- The "hessian" as is calculated is the derivative of the force imbalance. We can consider whether it is better to construct the hessian strictly as the second derivatives of the energy functional. Maybe the two slightly different numerical approaches will agree as the numerical resolution is increased.
- I have constructed a sequence of equilibria with "squeezed" current, and for each equilibrium I have calculated the eigenvalues of the hessian. When the eigenvalue becomes negative, there is a bifurcation! I think that this will be a really interesting case to explore, and we might soon have a paper describing how spec can determine stability. But, there is a lot I don't understand . . . I will upload some input files so that Joaquim can explore this.
- How should we name our input files? Soon, we will have many . . . For the test cases, I suggest something like G1V2Fi.001 for Lgeometry = 1 (G1), two volumes (V2), fixed boundary (Fi), and number 001. Perhaps we could add a date label. The input files that associated with publications can be placed in a suitably named subdirectory of InputFiles and can be named without any particular naming convention.
from spec.
Joaquim: I have uploaded two input files, G1V3Fi.xxx to master. These are two from a "squeezing" scan that I performed. The enclosed toroidal flux in the "squeezed" inner volume is constructed (using Echidna),
encflux = isequence * iscale
physicsprof[1].tflux = zero + ( one - encflux ) * half
physicsprof[2].tflux = one - ( one - encflux ) * half
physicsprof[3].tflux = one
Everything else in the input file is the same (except for output quantities such as the enclosed poloidal flux and helicity multiplier required to enforce the rotational-transform constraint and the geometry of the flux surfaces). This is very similar to the island squeezing that performed in measuring the 1/x and \delta function current densities. Do you remember how the island seemed to be of the opposite phase to what we expected?
I suggest you perform a similar scan, and for each case examine the minimum eigenvalue of the "hessian", as is calculated by spec in it's present form. The eigenvalue passes through zero as the toroidal flux becomes sufficiently squeezed, and the phase of the island flips. This must be associated with some type of instability, but I don't understand . . . When the current becomes too singular, does the equilibrium become unstable?
Please continue to study what I am calling the hessian, and please investigate this "instability" and how it is related to the negative eigenvalue. I think that, if we can understand this, we can easily write the first ever paper on linear stability for equilibria with islands, and this line of research will very much help us as we proceed onto the preconditioning, including the spectral constraints etc. in the general toroidal case.
from spec.
The "hessian" as is calculated is the derivative of the force imbalance. We can consider whether it is better to construct the hessian strictly as the second derivatives of the energy functional. Maybe the two slightly different numerical approaches will agree as the numerical resolution is increased.
OK, then why do you expect the Hessian matrix to be (tri)diagonal? Have a look at the last plots I attached: the Hessian is not (tri)diagonal.
How should we name our input files? Soon, we will have many . . . For the test cases, I suggest something like G1V2Fi.001 for Lgeometry = 1 (G1), two volumes (V2), fixed boundary (Fi), and number 001. Perhaps we could add a date label. The input files that associated with publications can be placed in a suitably named subdirectory of InputFiles and can be named without any particular naming convention.
Yes, I agree. But then let us create a "log-file", e.g. inputinfo.txt, that contains a one-sentence-explanation of what the file is (G,V,F are not all that counts, e.g. how is 001 different from 002, etc.). So for example each time someone adds an input file, e.g. G3V1Fi.001, then the person also modifies the file inputinfo.txt by adding one line with G3V1Fi.001 and some info.
What do you think?
from spec.
I have constructed a sequence of equilibria with "squeezed" current, and for each equilibrium I have calculated the eigenvalues of the hessian.
This is very similar to the island squeezing that performed in measuring the 1/x and \delta function current densities. Do you remember how the island seemed to be of the opposite phase to what we expected?
Please continue to study what I am calling the hessian, and please investigate this "instability" and how it is related to the negative eigenvalue. I think that, if we can understand this, we can easily write the first ever paper on linear stability for equilibria with islands, and this line of research will very much help us as we proceed onto the preconditioning, including the spectral constraints etc. in the general toroidal case.
@SRHudson, I agree that this is a very interesting study to carry out, and I am willing to do it, but I still need to convince myself that the hessian is not being messed up even in the 2-volume symmetric case (cartesian.sp you uploaded). I do not want to study a physical system with a tool that may have problems in exactly the objects (hessian) that we need in order to uncover the physics.
from spec.
I have been carrying out a verification of the Hessian in SPEC by comparing it to an exact analytical expression that one can get in a simple slab model. I wrote a two-pages pdf summarizing the results. Please have a look at it and let me know what do you think.
from spec.
I pushed the verification farther, by computing analytically the Hessian for a multi-volume slab equilibrium (with two continuous symmetries). A comparison with corresponding SPEC calculations indicate that SPEC is NOT computing the Hessian correctly!
Please have a look at the summary of this investigation:
from spec.
I have been looking in detail at dforce.h and it seems that by construction the matrix hessian(1:NGdof,1:NGdof)
should be tridiagonal in the slab case I considered. Looking at newton.h, however, where the hessian is written in the .DF file, I see that some manipulation on the hessian is performed,
ijdof = 0
do idof = 1, NGdof
do jdof = idof, NGdof ; ijdof = ijdof + 1 ; hessian(idof,jdof) = RR(ijdof) ! un-pack R matrix;
enddo
enddo
! derivative matrix = Q R;
hessian(1:NGdof,1:NGdof) = matmul( fjac(1:NGdof,1:NGdof), hessian(1:NGdof,1:NGdof) )
which I don't understand. @SRHudson What is conceptually being done here?
from spec.
Just a quick comment: please read through the documentation for hybrj. I think I had originally used the NAG routine C05PDF, and the documentation is online
https://www.nag.co.uk/numeric/FL/manual19/pdf/C05/c05pdf_fl19.pdf
In this documentation you will see "On final exit: the orthogonal matrix Q produced by the QR factorization of the final approximate Jacobian". So note that the Jacobian goes in as the Jacobian, but comes out as "Q". To get back to the Jacobian matrix, I had to multiply this by "R", which comes out packed as a vector. Note also that "The Jacobian is updated by the rank-1 method of Broyden." So, this routine may be changing the Jacobian. If you want the Jacobian, don't call this routine. Instead, first run SPEC to get to the solution. Then run SPEC to just calculate the Jacobian.
Does this make sense?
from spec.
The newton routine seems to be calling dforce.h to calculate the force given position and then calls the hybrj routine in order to find the zeros of the force. In doing so, it calculates some factorization matrices Q, R, that can allow retrieving the jacobian. And that is what seems to be done by newton.h right before the hessian is written on the .DF file.
What I do not understand is that the dforce.h subroutine itself calculates as well the hessian (and there I can explicitly see that by construction it should be tridiagonal), but does not seem to be used or at least does not seem to coincide with the written hessian in .DF
from spec.
I am not sure that I understand your question, but let me try to explain. Presently the PPPL computers have crashed and I cannot look at the source, so I am going by memory.
- The hessian is calculated in dforce.
- The hessian is used, and perhaps modified, in newton by hybrj. The potentially modified hessian is written to file. Most iterative routines for finding zeros will update the provided Jacobian using Broyden's method. If any iterations have been performed in newton/hybrj in order to find the solution, then it is likely that the Jacobian has been altered.
- A "clean" hessian is constucted by hessian.h as an aposteori diagnostic. There are options for comparing the hessian computed by dforce with a finite difference estimate and for computing the eigenvalues/vectors of the hessian, and these are related to the linear stability of the equilibrium. This would make a great paper, "The linear stability of partially relaxed MHD equilibria".
- What I am calling the hessian is really the derivative of the force vector, F = [[p+B^2/2]], with respect to geometry. I also sometime call this the derivative matrix.
I understand that the hessian as written to the .DF file by newton is slightly different from your analytic example, and my comments above suggest a reason why this might be the case. Can you remind me please: is SPEC getting the same solution as your analytic case for the geometry of the interfaces? I remember that we performed very careful investigations of SPEC vs analytic for our first paper on the singular currents, which was also in Cartesian. I can only guess that there is a misunderstanding at this point, or that a bug has been introduced recently. Either way, it should be easy to resolve this (perhaps by digging out older versions of SPEC if required).
Also, is there any chance that you (Joaquim) could spend a week or more here at PPPL? I am sure that we could resolve this matter quickly if we were working together (like the good old days).
from spec.
What I am calling the hessian is really the derivative of the force vector, F = [[p+B^2/2]], with respect to geometry. I also sometime call this the derivative matrix.
Yes, me too, and that is what I calculated also analytically.
If any iterations have been performed in newton/hybrj in order to find the solution, then it is likely that the Jacobian has been altered.
I see, but still the Hessian (as defined by your point 4) is unique - given position and constraints - and should not depend on what or how iterations are performed. I understand however that a numerical inversion to obtain the jacobian may pollute the retrieved hessian. If that is the case, then the .DF file does not contain the actual hessian and I should compare the analytical calculations to the hessian as calculated in dforce.h or in hessian.h. Right?
is SPEC getting the same solution as your analytic case for the geometry of the interfaces? I remember that we performed very careful investigations of SPEC vs analytic for our first paper on the singular currents, which was also in Cartesian.
As you can see in the document I uploaded (Hessian_verification.pdf), the geometry of the interfaces obtain in SPEC agrees well with theory but there is no knob I can turn that makes SPEC converge...in fact the only knob should be Lrad but this does change nothing! In fact using Lcheck=1, I can also see that the error in the Beltrami field does not change with Lrad...strange. In our first slab paper, we did not perform, technically speaking, a verification, in the sense that we did not show convergence of the error versus resolution. But it seemed to be the case that SPEC was retrieving the analytical results to high accuracy. What I could do is to perform this recent exercise with the old version of the code and see. I still find strange that Lrad does not change the error in beltrami field...
Also, is there any chance that you (Joaquim) could spend a week or more here at PPPL? I am sure that we could resolve this matter quickly if we were working together (like the good old days).
I like the good old days! I agree that we could work efficiently if I came ~1week. I will consult with the "supreme court" to see what I can do. If that is not an option in the short-term, we can try to find a week in which we both focus mostly on this problem and work together remotely. Let you know soon.
from spec.
Let's call SPEC1.0 the old version of SPEC that I have been running for all our papers so far, and that was verified to machine precision in our PoP2016.
Let's call SPEC2.0 the new version of SPEC (master branch) that we are working on.
As I reported in the last uploaded document above, a simple SPEC2.0 slab equilibrium calculation with only one inner interface and Mpol=Ntor=0 does not seem to retrieve the analytical results, not even as Lrad is increased (namely, neither the geometry of the interface nor the Hessian agree).
I tried running the same system with SPEC1.0 and found that as Lrad increases the geometry of the interface quickly converges to the theoretical estimate, easily approaching machine precision error.
This suggests that some sort of important bug has been introduced in SPEC2.0.
from spec.
Good.
I say good not because there is a bug but it is good that it has been identified that there is a bug!
Where can I get a copy of the source of SPEC1.0 ? It should be pretty easy to find the problem.
from spec.
It would be good to have SPEC1.0 uploaded on the GIThub together with (but decoupled from) the current SPEC source, some sort of tagged version of SPEC, named SPEC1.0. Is there a "best way" to do that?
Otherwise I can also send a zip file with the source, but it would be good to avoid old-school methods...Also, if we put SPEC1.0 on GIThub, then it would be pretty straighforward to reproduce all the verification calculations in slab,cylindrical,and stellarator geoemtries that we have carried out in recent papers.
from spec.
I have uploaded the source for SPEC1.0 by creating a folder within the master branch. In the future, we shall use "git tag" to tag stable versions of SPEC. Once the current version of SPEC is shown to be retrieving properly analytical results, shows convergence, the NAG branch is merged into master and tests performed, we can think of tagging SPEC as SPEC1.1, for example, with the idea that SPEC2.0 should be a verified free-boundary, non-stellarator symmetric, fast and robust code. Just an idea.
from spec.
Together with @SRHudson , we have identified that some bug was introduced in the code on Aug7th2017, when merging the ma00aa branch into master (the bug must then have been in that ma00aa branch). The log for this commit was:
commit 635915c
Merge: b98b8fe 5bf6a9d
Author: Joaquim [email protected]
Date: Mon Aug 7 14:20:38 2017 +0200
Merge pull request #14 from PrincetonUniversity/ma00aa
Confirming the merging of ma00aa branch into the master branch.
Checking out the master version before that commit, which can be done from the master branch by using:
git checkout b98b8fe
and then compiling and running this version, which produces, for example, the correct geometry of the interface in the slab 2-volume testcase I considered in my verification study (see document Hessian_verification.pdf uploaded below).
from spec.
Related Issues (20)
- GMRES error HOT 10
- Segmentation faults when compiled with cmake HOT 6
- Problem with master branch HOT 4
- Issue when restarting from .end file HOT 6
- Help needed to install python wrappers HOT 15
- Wrapper will not compile (reserved python words?) HOT 3
- Problem running tests on master branch HOT 3
- python_wrapper Makefile setup HOT 2
- Help needed for compilation with cmake HOT 8
- Question about force-gradient HOT 4
- Bug in force gradient when Lrzaxis=1 ?
- MATLAB metric subroutine needs fix HOT 2
- Python wrapper compilation for Henneberg representation branch HOT 2
- VMEC initializer HOT 2
- Read initial guess with python wrappers HOT 5
- Cannot install f90wrap HOT 7
- Question about matrix free method HOT 2
- SPEC license? HOT 5
- Installation on mac m1 HOT 16
- py_spec installation is borken HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
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.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from spec.