search

PrincetonUniversity / SPEC

The Stepped-Pressure Equilibrium Code, an advanced MRxMHD equilibrium solver.
https://princetonuniversity.github.io/SPEC/
GNU General Public License v3.0
25 stars 6 forks source link

Current constraint #75

Closed abaillod closed 4 years ago

abaillod commented 5 years ago

Hello everyone,

For the last months, I have been implementing a new constraint in SPEC - the current constraint. This will allow the user the specify the net toroidal current flowing in each volume and at each interface (as a current sheet).

This constraint is slightly different from the others because we don't evaluate it with quantities local to a volume, but global ones (more specifically the magnetic field on each side of the interface is needed to evaluate the current flowing at the interface). This required substantial changes in the file dforce.h to allow the hybrj routines to iterate on volumes solutions instead of one by one. If you are interested in more details, please refer to the document _Inputfiles/Verification/currentconstraint/TestCases_Comparison/SPEC_CurrentConstraintTestCaseVerification.pdf

I reached a point where I believe it is ready for merging with the master branch. If you have the time,

  1. Pull the CurrentConstraint branch and compile it. Then open Inputfiles/Verification/currentconstraint/TestCases_Comparison and follow the steps outlined in the ReadMe - this requires SPEC to run for about half an hour and compares the output of a local constraint with the current constraint (for the same equilibrium).

  2. Look at the plots - if the iota-profiles are different or the force difference between both runs is too large, report it here. Report as well any bugs / crashes / problems / questions here!

  3. If everything went smoothly, please notify me here as well. I won't merge with master until it has been tested on a few different computers.

Thank you!

PS: I write here some other minor changes that I implemented in the same branch

zhisong commented 5 years ago

@abaillod This is really great job! Well done! Is it convenient for you to talk about this work in detail in one of the SPEC meetings (potentially the next one)? We can go over your implementation in more detail and hopefully accelerate the merging. I notice there is conflicting files and we may need to work on resolving the conflict.

abaillod commented 5 years ago

Thank you Zhisong,

Yes talking about it during a SPEC meeting is a good idea - I will see with Joaquim when this could be possible.

I already tried (locally) a quick and dirty merge with master, to see how difficult it will be. It went rather smoothly - SPEC was working properly in fixed boundary. I couldn't test it in free boundary though...

abaillod commented 5 years ago

I merged the free_bound branch into my CurrentConstraint branch.

The current constraint (Lconstraint=3) is now working in free boundary as well. I created two input files for testing:

I noticed that in both cases the Lconstraint=3 calculation takes significantly more time to converge (more free boundary computations) than the Lconstraint=1 case. I still need to understand why...

jonathanschilling commented 4 years ago

@abaillod @jloizu I started evaluating the reference runs Antoine provided. A couple of these did not execute properly, i.e. took more than a few hours on the IPP Draco cluster or showed error messages, and I will re-execute them again one after another to figure out what exactly is going on.

abaillod commented 4 years ago

@jonathanschilling Thank you for taking the time to review this branch.

It is possible that some calculations with Lconstraint=3 take longer than the others. For now, only Lfindzero=1 is available with Lconstraint=3, i.e. the force gradient is evaluated numerically which is of course much slower than Lfindzero=2. However, test cases with Lconstraint=0,1,2 should converge in a matter of 10's of minutes.

jonathanschilling commented 4 years ago

@abaillod Thanks for the notice on the runtime. In fact, all runs I did so far completed in under a minute and I am curious to see which of the later ones stalls. I ran the test cases on the draco cluster at MPCDF with the following modules: intel/18.0.5, impi/2018.4, mkl/2018.4, git/2.16, hdf5-serial/1.10.5 and fftw-mpi/3.3.8

Here is my progress so far:

  1. G1V02L0Fi.001.sp, G1V02L0Fi.master.001.sp, G1V02L3Fi.001.sp These runs executed fine and the following output indicates that indeed the same equilibrium is computed:

    >> compare_spec_outputs('G1V02L0Fi.001.h5', 'G1V02L0Fi.master.001.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G1V02L3Fi.001.h5', 'G1V02L0Fi.master.001.h5');
Dabs:            5.5511e-16
Drel:            9.5651e-16
Estimate for df: 3.2216e-16
The two outputs can be considered the same

  2. G1V05L0Fi.001.sp, G1V05L0Fi.master.001.sp, G1V05L3Fi.001.sp When running G1V05L3Fi.001.sp, I got the following message: !!! WARNING: Using Lfindzero=2 with Lconstraint=3 might not converge. The hessian at fixed toroidal current has not yet been implemented. !!!. In your last comment, you mentioned that only for Lconstaint=3, Lfindzero=2 was not implemented yet? Apart from that, here is the output of the comparison routine:

    >> compare_spec_outputs('G1V05L0Fi.001.h5', 'G1V05L0Fi.master.001.h5');
Dabs:            3.3307e-16
Drel:            7.8127e-16
Estimate for df: 1.4199e-16
The two outputs can be considered the same
>> compare_spec_outputs('G1V05L3Fi.001.h5', 'G1V05L0Fi.master.001.h5');
Dabs:            1.0784e-12
Drel:            2.3061e-12
Estimate for df: 5.0424e-13
The two outputs are not the same

    Lconstraint=3 seems to deviate a little bit too much here; is this maybe related to the Lfindzero issue?

  3. G2V01L0Fi.001.sp, G2V01L0Fi.master.001.sp, G2V01L3Fi.001.sp In each of these runs, I get the following MKL error message:

    hesian :       0.78 : LHmatrix= T ;

Intel MKL ERROR: Parameter 4 was incorrect on entry to DGETRF.
hesian :            : myid=  0 ; idgetrf= -4 ; input error ; determinant=  1.00000E+00 ;

    For the run G2V01L3Fi.001.sp I get the warning related to Lfindzero=2 as above again.

    However, SPEC starts computing an equilibrium and indeed the outputs are the same:

    >> compare_spec_outputs('G2V01L0Fi.001.h5', 'G2V01L0Fi.master.001.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G2V01L3Fi.001.h5', 'G2V01L0Fi.master.001.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same

    Maybe this is an error and it is comparing zeros with zeros since the code crashed? I will check this and report back.

Also, the other comparisons are on the way.

abaillod commented 4 years ago

Hello Jonathan,

Answering point by point:

1) The warning !!! WARNING: Using Lfindzero=2 with Lconstraint=3 might not converge. The hessian at fixed toroidal current has not yet been implemented. !!! is normal. Using Lconstraint=3 with Lfindzero=2 will compute analytically the force gradient, using what has been implemented for the rotational transform constraint, i.e. Lconstraint=1.

For some reasons, it seems that it still converges towards the expected equilibrium faster than using Lfindzero=1; this is why I put a warning instead of throwing an error. This is of course temporary, and will be replaces as soon as I complete the implementation of the force gradient evaluation for Lconstraint=3.

2) Concerning the difference between G1V05L3Fi.001.sp and G1V05L3Fi.master.001.sp, this might be due to Lfindzero=2. I will try to run it with Lfindzero=1 to see if it reduces the difference.

3) About G2V01L0Fi.001.sp and similar, I think this is due to mistake in the input file. I checked on my computer and I get the same error - my bad! I will have a look.

Thank's a lot for the feedback!

abaillod commented 4 years ago

I solved point 3. - it was indeed due to a mistake in the input file. SPEC was attempting to evaluate the hessian, though NGdof was set to 0. I changed LHmatrix from T to F in the input files.

Maybe this should be tested in general at the beginning of hessian.f90... I don't know if it makes sense to call this subroutine with NGdof=0.

abaillod commented 4 years ago

@jonathanschilling, I looked into point (2).

You were right, using Lfindzero=2 with Lconstraint=3 was leading to a slightly different equilibrium. Lfindzero=2 lead to a difference in force balance of df\~1E-12 in comparison to df\~1E-14 for Lfindzero=1.

I thus changed preset.f90 in order to raise an error if Lconstraint=3 is used with Lfindzero=2. This is of course temporary until I figure out how to implement the force gradient (hopefully in February).

I modified the input files in the Verification directory accordingly: all L3 input file now use Lfindzero=1.

Please pull the new version and tell me if everything works out for you! I think I addressed all issues you raised in your last comment.

jonathanschilling commented 4 years ago

@abaillod Thanks for adressing the issues. I will re-run the input files and report back about the results.

jonathanschilling commented 4 years ago

@abaillod I finished running all test cases provided by you today and would like to report on the results. Briefly speaking, all but the last test case performed very well now and I am positive that the merge is approaching fast :-) Thanks for your improvements!

  1. The case G1V02L0Fi.001 still runs fine and produces the same deviations as before.

    >> compare_spec_outputs('G1V02L0Fi.001.h5', 'G1V02L0Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G1V02L3Fi.001.h5', 'G1V02L0Fi.master.001.sp.h5');
Dabs:            5.5511e-16
Drel:            9.5651e-16
Estimate for df: 3.2216e-16
The two outputs can be considered the same

  2. The case G1V05L0Fi.001 does not produce any warnings anymore and the deviations are much reduced:

    >> compare_spec_outputs('G1V05L0Fi.001.h5', 'G1V05L0Fi.master.001.sp.h5');
Dabs:            3.3307e-16
Drel:            7.8127e-16
Estimate for df: 1.4199e-16
The two outputs can be considered the same
>> compare_spec_outputs('G1V05L3Fi.001.h5', 'G1V05L0Fi.master.001.sp.h5');
Dabs:            1.5155e-14
Drel:            3.5548e-14
Estimate for df: 6.4606e-15
The two outputs can be considered the same

  3. The case G2V01L0Fi.001 also does not produce any warnings/error messages anymore and the deviations are fine:

    >> compare_spec_outputs('G2V01L0Fi.001.h5', 'G2V01L0Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G2V01L3Fi.001.h5', 'G2V01L0Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same

  4. The case G3V01L0Fi.001 ran fine and the deviations are ok:

    >> compare_spec_outputs('G3V01L0Fi.001.h5', 'G3V01L0Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G3V01L3Fi.001.h5', 'G3V01L0Fi.master.001.sp.h5');
compare_spec_outputs('G3V01L3Fi.001.h5', 'G3V01L0Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same

  5. The case G3V01L0Fi.002 ran fine and the deviations are ok:

    >> compare_spec_outputs('G3V01L0Fi.002.h5', 'G3V01L0Fi.master.002.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G3V01L3Fi.002.h5', 'G3V01L0Fi.master.002.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same

  6. The case G3V02L1Fi.001 ran fine and the deviations are ok:

    >> compare_spec_outputs('G3V02L1Fi.001.h5', 'G3V02L1Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G3V02L3Fi.001.h5', 'G3V02L1Fi.master.001.sp.h5');
Dabs:            8.5265e-14
Drel:            9.3566e-09
Estimate for df: 7.77e-19
The two outputs can be considered the same

  7. The case G3V08L1Fi.001 stalled for more than 30 minutes at the calculation of the force gradient matrix when run with 1 CPU. I increased the number of CPUs to 8 and kept it like this for all following runs. Then the case executed fine and the deviations are ok:

    >> compare_spec_outputs('G3V08L1Fi.001.h5', 'G3V08L1Fi.master.001.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G3V08L3Fi.001.h5', 'G3V08L1Fi.master.001.sp.h5');
Dabs:            6.0432e-12
Drel:            9.8166e-06
Estimate for df: 3.7202e-18
The two outputs can be considered the same

  8. The case G3V08L1Fi.002 executed fine with 8 CPUs and the deviations are ok:

    >> compare_spec_outputs('G3V08L1Fi.002.h5', 'G3V08L1Fi.master.002.sp.h5');
Dabs:            0
Drel:            0
Estimate for df: 0
The two outputs are exaclty the same
>> compare_spec_outputs('G3V08L3Fi.002.h5', 'G3V08L1Fi.master.002.sp.h5');
Dabs:            7.514e-13
Drel:            9.2887e-05
Estimate for df: 6.0784e-21
The two outputs can be considered the same

  9. The case G3V08L1Fr.001 executed fine with 8 CPUs for the master branch and the local constraint. The Lconstraint=3 took much longer to execute (+1h) and the deviations are significant:

    >> compare_spec_outputs('G3V08L1Fr.001.h5', 'G3V08L1Fr.master.001.sp.h5');
Dabs:            1.45e-13
Drel:            5.4589e-09
Estimate for df: 3.8512e-18
The two outputs can be considered the same
>> compare_spec_outputs('G3V08L3Fr.001.h5', 'G3V08L1Fr.master.001.sp.h5');
Dabs:            0.00070288
Drel:            0.48995
Estimate for df: 1.0084e-06
The two outputs are not the same

Here is my logfile so far: log_current_constraint_review.txt

I am actually a little bit irritated that the deviations are zero for so many cases. I assumed that due to the different constraints we would get deviations on the order of numerical accuracy. @abaillod , @jloizu , @SRHudson : What do you think about this? Do we have to worry that the tests are not sufficient yet to determine the (expected?) deviations between the runs?

abaillod commented 4 years ago

@jonathanschilling

Thank you for all these tests!

It would be worth it to check if the case G3V08L1Fi.001 also stalls when run with the master branch - to be sure that my branch doesn't impact too much SPEC performances. I don't really mind if Lconstraint=3 is slow; however Lconstraint=1 should be as fast as before... I will also test this on my local computer.

Concerning the cases where |df|=0: this is normal when only one volume is considered - no force balance is required on the interfaces (since there are no interfaces). I actually use the matlab routine _compare_specoutputs like a black box - I hope the test is sufficient to say if two equilibria are equivalent or not.

Finally the differences in G3V08L1Fr.001: I think this is due to the tolerance in Bn - I am running both cases again with a lower tolerance, I will tell you if it improves or not.

Thank's again for your help!

abaillod commented 4 years ago

To answer to my last comment:

It would be worth it to check if the case G3V08L1Fi.001 also stalls when run with the master branch - to be sure that my branch doesn't impact too much SPEC performances. I don't really mind if Lconstraint=3 is slow; however Lconstraint=1 should be as fast as before... I will also test this on my local computer.

I ran this case with the current constraint branch and the master branch. Both took approximately the same amount of time to complete - with one or more CPUs.

Concerning the cases where |df|=0: this is normal when only one volume is considered - no force balance is required on the interfaces (since there are no interfaces). I actually use the matlab routine _compare_specoutputs like a black box - I hope the test is sufficient to say if two equilibria are equivalent or not.

To test if the field in a volume is the same, I implemented an additional matlab routine called _compare_specoutputs2.m. This routine checks that the difference in each Fourier harmonic of the field and geometry between two files is below a given tolerance. I tested all V01 cases and all passed the test, excepted the free boundary one (see below)

Finally the differences in G3V08L1Fr.001: I think this is due to the tolerance in Bn - I am running both cases again with a lower tolerance, I will tell you if it improves or not.

This is not due to Bntol - decreasing it does not reduce the error. Actually, the comparison is non-trivial: it is rather hard to converge to the same equilibrium given different constraints! To show that, I compared two equilibria obtained once with Lconstraint=0 and once with Lconstraint=1 using the master branch. They didn't pass the test. It basically means that free boundary equilibria obtained with different constraints won't be the same up to machine precision - thus G3V08L1Fr and G3V08L3Fr are not comparable. However, note that both the master and the current constraint branch successfully find the same equilibrium for G3V08L1Fr, so the current constraint branch has no effect on the free boundary implementation for previously existing constraints (L0 and L1).

With these points addressed, I think all test cases passed successfully at IPP and SPC? Please remind me if I forgot a bug / problem!

smiet commented 4 years ago

Since my latest commit to master, you will see that the file numrec.f90 is conflicting. This is because the file had non-unix line endings, which I have reverted. Please overwrite this file; git checkout master -- numrec.f90 should do the trick. (or merge master into your branch; should give no other conflicts.)

(this was caught by the linting step in the newly set-up regression tests. I will be adding physics tests, but for now it can catch simple problems like this. Merge master into your branch to see it in action!)

abaillod commented 4 years ago

Thank you @smiet for noticing the conflict, I corrected it following your indication. I tested again the file G3V08L3Fi.001.sp and everything works as before.

zhisong commented 4 years ago

I came up with an idea. Again, this is just a suggestion. You can choose to implement it if want, but it is not necessary given that your code is working now.

In fact, if all you need to iterate is pflux and the only thing you care about is Bt00 the "averaged" toroidal field on the surface, then there is a much simpler way to get the current constraint without any iteration (hybrd1 or something).

Bt00 is a linear combination in the field vector, Atemn, Azemn. What I mean by linear combination is that, if we let a=(Atemn, Azemn,some Lagrange multipliers)^T, a total of NAdof degrees of freedom, we will get Bt00=M dot a, in which M is some matrix with size 1*NAdof (a row vector).

On the other hand, your field vector is solved by a=(dMA- mu * dMD)^-1 dot dMB dot (psi_t, psi_p)^T. Considering above, we will get Bt00=M dot (dMA- mu * dMD)^-1 dot dMB dot (dtflux, dpflux)^T = alpha dtflux + beta dpflux, where alpha and beta are two scalar constants. That means, the relationship between Bt00 and dtflux and dpflux is linear.

Using this relationship, you can cast the current constraint as a set of linear equations with unknows dpflux in each volume. This can be solved once after the field in all the volumes are computed. No iteration is needed. And you don't need to store the dMA matrix since all you need is the coefficients alpha and beta. After you found the right dpflux that satisfying the current constraint, the solution a can be reassembled by a linear combination without going through the solving process again. I would expect the whole process takes approximately the same amount of time as the Lconstraint=0 case.

To compute the force gradient, you will need to differentiate the row vector M above with respect to the interface geometry.

abaillod commented 4 years ago

@Zhisong: Thank you for the idea! Yes, that might speed up the code but it would require to re-implement completely the constraint. I think it is better to first merge this branch once it is working and in a second time, when we have the time, we can try to optimize it. This is definitely something to keep in mind in case we notice that the constraint is very slow.

Speaking of merge, I think all detected problems have been solved (automatic tests have passed). @zhisong, could you try to run the tests on your computer if you didn't do it already? Try to compile the branch and run it on each case in the ci/ directory (that should be enough).

smiet commented 4 years ago

Speaking of a merge, let's add a CI test that put's the new features through it's paces on every new update to Master!

@abaillod , in what geometry shall we test the currentconstraint? Could you provide one or several .sp input file(s) that puts the new features through it's paces, as well as a .h5 output file of the run that you have checked gives physically reasonable answers?

I will then run the cases in the docker image, compare to your h5, and add the cases to the CI so all future merges to master keep the CurrentConstraint functionality.

abaillod commented 4 years ago

Good point @smiet! I created two folders in /ci: one for the case G3V08L3Fi and one for the same case in free boundary, G3V08L3Fr. Both have the input file and a compare.h5 file. They both have a reasonable resolution such that the execution time is small (~100s on a single CPU).

I am not sure how to modify the docker file though, I let you do the required modifications!

jloizu commented 4 years ago

I think that after @zhisong tests and approves, we should be able to merge, finally, this branch into master.

Please @jonathanschilling and @zhucaoxiang approve merge if you think it is appropriate.

smiet commented 4 years ago

I have added to the CI and it should now automatically test the two CurrentConstraint tests. When reproducing on my own laptop, the fixed boundary ran fine, but I got a mismatch in the freeboundary solution.

the unmatched variables are: UNMATCHED: volume , diff= 3.90799E-12

Elements in poincare UNMATCHED: R , diff= 3.64714E-09 UNMATCHED: Z , diff= 5.96654E-09 ok: s element average difference = 2.4220016305655196e-14 ok: success element average difference = 0.0

Elements in transform ok: diotadxup element average difference = 0.0 UNMATCHED: fiota , diff= 4.58015E-11

full log file: CurrentConstraint_free_log.txt

Now that the Poincar\'e is different is unsurprising (and I will remove that from the regression testing as even minute differences blow up when integrating in a chaotic field for a long distance), but volume and fiota are a little more worrying. Though the differences are very small, and can maybe be ignored.

@jloizu, @abaillod, do you think a difference of this magnitude is small enough to be irrelevant? In that case I will replace the compare.h5 file with one generated by docker, so the tests will pass as they are generated almost identically.

Then we should be good to merge!

zhisong commented 4 years ago

I tried all the cases in the ci folder and ran the tests by hand. I ran it on two machines: my Mac with gfortran and the local cluster GADI with ifort. Here is what I get.

In the existing cases, all fine except on the cluster, the beltramierror does not match. On my Mac, there is no such problem.

I think this is fine, since the Beltrami field error depends on the local LAPACK library and is machine specific. I would suggest @smiet to remove the comparison for beltramierror. Maybe it is fine to unset Lcheck=1 in the input files?

In the new cases:

G3V08L3Fi Cluster: UNMATCHED: IPDt , diff= 1.38889E-02 Mac: UNMATCHED: IPDt , diff= 1.11111E-01

G3V08L3Fr Same as @smiet . They are at different values on the cluster and on my Mac. I think this may be a consequence of virtual casing. @smiet is it possible to add a case with virtual casing to the list of ci? Like some inputs in the VMEC-SPEC comparison.

abaillod commented 4 years ago

I compiled SPEC with gcc to run all the tests a second time (I didn't do that before pushing the branch). Results are:

G3V08L3Fi: Everything matches

G3V08L3Fr: 

UNMATCHED: R , diff= 5.90969E-09
UNMATCHED: Z , diff= 8.28180E-09
UNMATCHED: fiota , diff= 1.03251E-10

I also think that the Poincare integration is not a big deal. However I don't get why fiota is different when computed with intel instead of gcc. Is this related to the Poincare integration as well? @jloizu maybe you have an idea?

@zhisong It is very weird that you get such a large difference in IPDt. Could you try to run the test G3V08L3Fi again with the latest version of the branch, and send me the log? And how do you compile SPEC?

jloizu commented 4 years ago

I compiled SPEC with gcc to run all the tests a second time (I didn't do that before pushing the branch). Results are:

G3V08L3Fi: Everything matches

G3V08L3Fr:

UNMATCHED: R , diff= 5.90969E-09
UNMATCHED: Z , diff= 8.28180E-09
UNMATCHED: fiota , diff= 1.03251E-10

I also think that the Poincare integration is not a big deal. However I don't get why fiota is different when computed with intel instead of gcc. Is this related to the Poincare integration as well? @jloizu maybe you have an idea?

Well yes the calculation of iota inside the volumes in SPEC is done via the field-line-tracing which is carried out after the run in order to produce Poincare trajectories (and iota). So I can imagine that indeed if R,Z from the poincare values are different, that iota might be as well.

jonathanschilling commented 4 years ago

I am currently reviewing. For this, I execute again all the runs in TestCases/Verification/CurrentConstraint with the latest state of the code, both on 1 CPU and on 32 CPUs, both using Intel and GCC compilers. Once they are complete, I will let you know about the results.

zhisong commented 4 years ago

All clear for me.

jloizu commented 4 years ago

@jonathanschilling please provide some feedback on this pull request

@abaillod I suggest that in a couple of days you just merge the branch into master in any case

abaillod commented 4 years ago

To summarize:

I think we thus all agree and I can (finally) proceed with this merge! Thank you to all who participated and tested the branch. I will now drink (a) beer(s) to celebrate...

zhucaoxiang commented 4 years ago

@abaillod I think we may need a quick fix to resolve some minor issues in the master. See the following snippet, it contains unresolved merge conflicts.

https://github.com/PrincetonUniversity/SPEC/blob/master/Utilities/pythontools/compare_spec_BACKUP_5088.py#L8-L12

abaillod commented 4 years ago

@zhucaoxiang True, thank's for spotting the error.

I am not using the python tools and for some reasons I can't import the module when running python, so that makes it difficult for me to debug and test that. What is the correct way to import py_spec, from py_spec.read_spec import SPEC or from py_spec import SPEC?

Anyway I will open soon a new pull request for the ForceGradient_CurrentConstraint branch - we can discuss and fix that with this future merge.