Failed to run example case 1D-neon-source

[lqs661]

When I tried to run the 1D-neon-source case in example/tokamak-1D/extra, I met some problems. The log indicated the following error:
#0  (filtered)
#1  (filtered)
#2  in Options& set(Options&, Field3D) [clone .isra.0] [clone .cold]
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/hermes-3
#3  in decltype(auto) add<Field3D, GuardedOptions, void>(GuardedOptions&&, Field3D)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/hermes-3
#4  in NeutralParallelDiffusion::transform_impl(GuardedOptions&)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/libhermes-3-lib.so
#5  in Component::transform(Options&)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/libhermes-3-lib.so
#6  in ComponentScheduler::transform(Options&)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/libhermes-3-lib.so
#7  in Hermes::rhs(double, bool)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/hermes-3
#8  in Solver::run_rhs(double, bool)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/external/BOUT-dev/lib/libbout++.so.5.2.0
#9  in SNESSolver::rhs_function(_p_Vec*, _p_Vec*, bool)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/external/BOUT-dev/lib/libbout++.so.5.2.0
#10 in SNESSolver::snes_function(_p_Vec*, _p_Vec*, bool)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/external/BOUT-dev/lib/libbout++.so.5.2.0
#11 in SNESComputeFunction
    at /public3/home/t6s010954/library_mine/petsc/src/snes/interface/snes.c:2496:7
    2494:       SNESFunctionFn *computefunction;
    2495:       PetscCall(DMSNESGetFunction(dm, &computefunction, &ctx));
  > 2496:       PetscCallBack("SNES callback function", (*computefunction)(snes, x, y, ctx));
                ^
    2497:     }
    2498:     PetscCall(VecLockReadPop(x));
#12 in SNESLineSearchApply_BT
    at /public3/home/t6s010954/library_mine/petsc/src/snes/linesearch/impls/bt/linesearchbt.c:141:7
     139:       PetscCall(SNESComputeObjective(snes, W, &g));
     140:     } else {
   > 141:       PetscCall((*linesearch->ops->snesfunc)(snes, W, G));
                ^
     142:       if (linesearch->ops->vinorm) {
     143:         gnorm = fnorm;
#13 in SNESLineSearchApply
    at /public3/home/t6s010954/library_mine/petsc/src/snes/linesearch/interface/linesearch.c:643:3
     641:   PetscCall(PetscLogEventBegin(SNESLINESEARCH_Apply, linesearch, X, F, Y));
     642:
   > 643:   PetscUseTypeMethod(linesearch, apply);
            ^
     644:
     645:   PetscCall(PetscLogEventEnd(SNESLINESEARCH_Apply, linesearch, X, F, Y));
#14 in SNESSolve_NEWTONLS
    at /public3/home/t6s010954/library_mine/petsc/src/snes/impls/ls/ls.c:234:5
     232:        and evaluate F = function(X) (depends on the line search).
     233:     */
   > 234:     PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, Y));
              ^
     235:     PetscCall(SNESLineSearchGetReason(linesearch, &lssucceed));
     236:     PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
#15 in SNESSolve
    at /public3/home/t6s010954/library_mine/petsc/src/snes/interface/snes.c:4905:5
    4903:     snes->reason = SNES_CONVERGED_ITERATING;
    4904:     PetscCall(PetscLogEventBegin(SNES_Solve, snes, 0, 0, 0));
  > 4905:     PetscUseTypeMethod(snes, solve);
              ^
    4906:     PetscCall(PetscLogEventEnd(SNES_Solve, snes, 0, 0, 0));
    4907:     PetscCheck(snes->reason, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Internal error, solver %s returned without setting converged reason", ((PetscObject)snes)->type_name);
#16 in SNESSolver::run()
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/external/BOUT-dev/lib/libbout++.so.5.2.0
#17 in Solver::solve(int, double)
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/external/BOUT-dev/lib/libbout++.so.5.2.0
#18 in main
    at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/hermes-3
#19 (filtered)
#20 at /public3/home/t6s010954/hermes-dev/hermes-3-boron-test/build/hermes-3

====== Exception thrown ======
Setting invalid value for 'species:ne:density_source'

This is the first time I have encountered this error. I have tried to change the initial values, remove the source of the neutrals of Neon, and even reduce Neon’s charge state to +4, yet none of these measures proved effective. Could you please explain what might be causing it? Thank you.

[mikekryjak]

Sorry, we haven't looked at some examples for a while and they need checking. This error message probably means there is a NaN in the electron density source. It could be due to some division by zero error due to some part of the system not being floored properly. We have several PRs open on improving the floors further.

This is likely happening near the start of the simulation as it's initialised with flat profiles and so there will be many shocks. I would try to start from an example without impurities as an initial condition.

[lqs661]

I have tried to start the case from a pure Deuterium case as the initial condition. Besides, I simplified the settings, which just included only neutral of neon and ne+1. The option 'timestep' has also been reduced to 500. The source of 'Nne' has been removed. However, the case still failed to converge. The problem were the same "Setting invalid value for 'species:ne:density_source'". Is this program not suitable for solving the complete one-dimensional transport equation of impurities?

[mikekryjak]

@lqs661 Hermes-3 is a research code in development, and there haven't yet been serious attempts to perform impurity transport simulations, although the functionality for this is complete - apart from one group that has just started looking at carbon. We are implementing better parallel closure this year and will be attempting argon simulations, but this has not yet begun. It's on my list of things to do to review all of the examples, but I need to do several other things first I'm afraid. Have you checked for any NaNs in the solution?

[bendudson]

Hi @lqs661 ! I'm looking into this and will update the examples. The code changes over time and unfortunately we do not currently have an automated way to test that the examples still work.

Please try setting the diffusion collisions mode to "multispecies":

[neutral_parallel_diffusion]
diffusion_collisions_mode = multispecies

The default is "afn" that only includes ionization and charge-exchange. This doesn't include neutral-neutral collisions, so when the plasma density goes to zero the collision frequency (nu) goes to zero and so Dn is NaN. Setting diffusion_collisions_mode to multispecies restores the original collision frequency calculation that includes all collisions.

[bendudson]

Hi @lqs661 Please see PR #553 and try the examples in the fix-1d-examples branch.

[lqs661]

Well, I have tried the fixed 1D-neon case. The running still failed, and the error messages as follows:
======== SNES failed =========

Return code: 92, reason: 0
Nd+ : (1.0000002160704649 -> 6.268656619975728), ddt: (-1.0439702703422534e-05 -> 1.2838966876593805e-05)
Pd+ : (0.9696900801053795 -> 1.77773510720024), ddt: (-0.007719844380358699 -> 0.0013483174429669439)
NVd+ : (-1.0666680170036139e-07 -> 3.6754342116671226), ddt: (-0.07189996430957274 -> 1.2226269551647026e-06)
Nd : (2.4628467311274275e-07 -> 4.391306619342075), ddt: (-0.04423960278295617 -> 0.16589652884081307)
Pd : (3.670344816286899e-07 -> 0.5480353167831461), ddt: (-0.0018184919765497258 -> 0.008783389112977415)
NVd : (-0.014995735820001656 -> 0.1925953962754313), ddt: (-0.0012078342278194591 -> 0.02783805218866205)
Nne : (0.01 -> 0.01), ddt: (-5.079673791284105e-05 -> -3.249677556667195e-05)
Pne : (0.0001 -> 0.0001), ddt: (-3.7937195216784395e-07 -> 2.5016775459168515e-07)
NVne : (0 -> 0), ddt: (-1.4913339393748115e-07 -> 1.3240927275237145e-06)
Nne+ : (1e-05 -> 1e-05), ddt: (3.285563651819695e-05 -> 5.079110604472144e-05)
Pne+ : (1e-07 -> 1e-07), ddt: (3.3004826338742405e-07 -> 5.093865893380955e-07)
NVne+ : (0 -> 0), ddt: (-2.9611387368856915e-08 -> 3.204688261060151e-07)
Nne+2 : (1e-05 -> 1e-05), ddt: (-1.8417255250129475e-08 -> 2.460955090980972e-08)
Pne+2 : (1e-07 -> 1e-07), ddt: (7.397506840418423e-10 -> 6.800808689910927e-09)
NVne+2 : (0 -> 0), ddt: (-1.262095366006504e-07 -> 1.4316434528940686e-07)
Nne+3 : (1e-05 -> 1e-05), ddt: (-2.6160645944879197e-08 -> 1.0141239644671576e-08)
Pne+3 : (1e-07 -> 1e-07), ddt: (1.1231130247720677e-09 -> 1.439825956312947e-08)
NVne+3 : (0 -> 0), ddt: (-2.951439343407759e-07 -> 7.236896319099873e-08)
Nne+4 : (1e-05 -> 1e-05), ddt: (-2.6256207723263735e-08 -> 9.925215823563905e-09)
Pne+4 : (1e-07 -> 1e-07), ddt: (1.907328967776614e-09 -> 2.48257963935076e-08)
NVne+4 : (0 -> 1.0916187057447417e-07), ddt: (-5.355857227857426e-07 -> -5.0915860733111344e-11)
Nne+5 : (1e-05 -> 1e-05), ddt: (-2.613580386073639e-08 -> 1.0000134905189356e-08)
Pne+5 : (1e-07 -> 1e-07), ddt: (2.918472532319436e-09 -> 3.789863814121506e-08)
NVne+5 : (0 -> 0), ddt: (-8.479956193209592e-07 -> -8.03257655120191e-11)
Nne+6 : (1e-05 -> 1e-05), ddt: (-2.602158264244017e-08 -> 1.0080990076838554e-08)
Pne+6 : (1e-07 -> 1e-07), ddt: (4.1375625481002944e-09 -> 5.354618709841034e-08)
NVne+6 : (0 -> 0), ddt: (-1.2328593566751348e-06 -> -1.1640132284262145e-10)
Nne+7 : (1e-05 -> 1e-05), ddt: (-2.5886348476508823e-08 -> 1.0194739504002126e-08)
Pne+7 : (1e-07 -> 1e-07), ddt: (5.5633158383542405e-09 -> 7.171339827428833e-08)
NVne+7 : (0 -> 0), ddt: (-1.6905263968227222e-06 -> -1.5914253396572911e-10)
Nne+8 : (1e-05 -> 1e-05), ddt: (-2.5734006907717418e-08 -> 1.0336380376784008e-08)
Pne+8 : (1e-07 -> 1e-07), ddt: (7.188349501530235e-09 -> 9.235340744308864e-08)
NVne+8 : (0 -> 0), ddt: (-2.221292553195217e-06 -> -2.08549399995845e-10)
Nne+9 : (1e-05 -> 1e-05), ddt: (-2.563154615609484e-08 -> 1.0421226620848042e-08)
Pne+9 : (1e-07 -> 1e-07), ddt: (9.009319611435154e-09 -> 1.1542473276213168e-07)
NVne+9 : (0 -> 0), ddt: (-2.8254146478715693e-06 -> -2.646219219593447e-10)
Nne+10 : (1e-05 -> 1e-05), ddt: (-1.5861755456323852e-07 -> -5.411708996377302e-12)
Pne+10 : (1e-07 -> 1e-07), ddt: (1.1023206675105308e-08 -> 1.3981225381739696e-07)
NVne+10 : (0 -> 0), ddt: (-3.5031199022907283e-06 -> -3.2736010081722774e-10)
Pe : (0.6628840587563763 -> 0.9324370199608221), ddt: (-0.00022546587715885358 -> 0.0002236946930133793)

Besides, I have tried another method to running the neon case. First, I just remained the ne+0 to ne+4, removed the 'braginskii_ion_viscosity' and set the dneut=0.01. Restarting from the initial state of pure D case, the 'reduced' neon case has run successfully. Then I add the higher charge species gradually until to ne+10. Third, I increased the dneut from 0.01 to 10 gradually, setting it at 0.01, 0.05, 0.1, 0.5, 1.0, 5.0, 10.0. The case could also run successfully. Of course, the whole process span a long time. However, I have not tried to add the 'braginskii_ion_viscosity'.

[lqs661]

Recently, I have tried the final test of 1D-neon case. The test results were less than satisfactory. I added the 'braginskii_ion_viscosity' component and switch on 'thermal_conduction' for all of neon particles, including ne+0 to ne+10. The input file and results are as follows. As the picture shown, the case failed to converge. The target temperature oscillates over time. I tried to increase the timestep to ensure the convergence of the case, but it was unsuccessful.

BOUT_inp.txt

[mikekryjak]

Hi @lqs661. A few comments/questions:

• Are you running with the latest master branch?
• The timestep setting is just the user output frequency, and has no impact on convergence. The real solver timestep is adaptive and controlled by the solver.
• It looks like while there is a oscillation, the case profiles are still evolving. Have you run it for longer? How long is it taking to run?
• I suspect that your issues could be due to very low densities or pressures of some of the neon species. The first thing to try is to see if you are hitting the density or temperature floor, and if so then reduce them:

  hermes-3/src/evolve_pressure.cxx

  Lines 42 to 51 in 44ce2e6

  	density_floor = options["density_floor"].doc("Minimum density floor").withDefault(1e-7);
  	low_n_diffuse_perp = options["low_n_diffuse_perp"]
  	.doc("Perpendicular diffusion at low density")
  	.withDefault<bool>(false);
  	temperature_floor = options["temperature_floor"]
  	.doc("Low temperature scale for low_T_diffuse_perp")
  	.withDefault<BoutReal>(0.1)
  	/ get<BoutReal>(alloptions["units"]["eV"]);

  . Note that the floor units are normalised.
• Alternatively, you can try the new more consistent flooring scheme (Fluid equation of state with density floor #359) to see if it eliminates the fluctuations. But ideally you don't want to cross the floor.