robust convergence

[jloizu]

Following here the discussion that started with the iter-free-boundary testcase, which does not "robustly converge".

Below, I copy a comment by @SRHudson that summarizes this problem:

I think that the failure to robustly converge for this free-boundary iter-like case is associated with regularization problem near the coordinate/magnetic axis. The geometry of innermost interface goes crazy. Please see http://w3.pppl.gov/~shudson/Spec/packxi.pdf where I have given a brief description of the coordinate pre-conditioning factor that I have included. Also, I think that the spectral constraints, http://w3.pppl.gov/~shudson/Spec/lforce.pdf are an absolute mess. The spectral constraints need to be revised. I think I know what to do, and I think that I am close to having this problem fixed once and for all, but I have not had the opportunity to concentrate on this problem. If this matter is deemed a priority, then I will work on this. Note that I can get spec to converge for this case. It just requires a little human intervention and playing around with different choices of Linitialize. See http://w3.pppl.gov/~shudson/Spec/global.pdf for the possible values of Linitialize. There is, however, potentially a more serious problem with this iter-like free-boundary case. Do I have the correct sign on the plasma current? What is the rotational transform negative? Do I have the correct normalization on the magnetic field at the computational boundary produced by the plasma current (computed using the virtual casing)? I uploaded this input file for reference. It would be a lot easier to verify that free-boundary SPEC is working for tokamaks using a circular cross section boundary configuration.

[jloizu]

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.

[jloizu]

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

[SRHudson]

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

[jloizu]

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?

[SRHudson]

1) I pulled the master and found the input files, thanks. Are the input files available on other branches? 2) 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.)

[SRHudson]

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

[jloizu]

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.

iterfixed_poincare_vessel.pdf

[jloizu]

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;

[SRHudson]

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.

[SRHudson]

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.

[jloizu]

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.

[jloizu]

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?

[jloizu]

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?

[SRHudson]

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!

[jloizu]

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.

[jloizu]

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.

[SRHudson]

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

[jloizu]

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?

[jloizu]

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?

[SRHudson]

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.

[SRHudson]

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.

[jloizu]

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?

[SRHudson]

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.

[jloizu]

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

circulartok_poincare.pdf

[jloizu]

(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?

[SRHudson]

Joaquim: 1) 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. 2) 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. 3) 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.) 4) 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. 5) 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.

[SRHudson]

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

[jloizu]

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.

[zhucaoxiang]

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.

[SRHudson]

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.

[jloizu]

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.

[jloizu]

For the circular tokamak equilibrium, here is what I get for the Hessian matrix.

circulartok_hessian1.pdf

circulartok_hessian2.pdf

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.

[jloizu]

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!

[zhucaoxiang]

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.

[SRHudson]

1) I suggest you simplify the equilibrium as much as possible. Perhaps you can use only one internal interface and low Fourier resolution. 2) 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. 3) 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. 4) 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. 5) Start with Cartesian. There are no spectral constraints, the geometric deformations are easy to identify, . . . 6) Just investigate everything that you can!

[SRHudson]

1) 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. 2) I have added an input file in Cartesian geometry. 3) 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. 4) 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). 5) 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. 6) 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? 7) 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. 8) 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. 9) 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. 10) 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).

[jloizu]

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.

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

[jloizu]

I fixed this problem by doing

git checkout 075c256fbc9f7e4d 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! :)

[jloizu]

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)

[jloizu]

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,

cartesian_hessian.pdf

(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?

[jloizu]

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?

[SRHudson]

"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

[jloizu]

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.

[jloizu]

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.

[jloizu]

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.

[SRHudson]

1) 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. 2) 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. 3) 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.

[SRHudson]

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.

[jloizu]

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?

[jloizu]

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.

[jloizu]

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.

Hessian_verification.pdf