ERROR KGslErrorHandler.cxx and domain error (1) in ellint.c

[marcogobbo]

I'm trying to simulate an electron gun with a simple setup that includes two electrodes with one with a hole. The electrons generated from the cathode are accelerated to the anode with a hole. Two termination events are defined ksterm_max_z and ksterm_death, this one seems to produce this error each time the electrons collide with the anode surface

16:55:04 [ERROR] l/Utility/KGslErrorHandler.cxx:84    GSL error: ellint.c:542: domain error (1)
  1# gsl_sf_ellint_Ecomp_e
  2# gsl_sf_ellint_Ecomp
  3# KEMField::KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing(double const*, double const*)
  4# KEMField::KGaussianQuadrature::operator()(double (**)(double const*, double const*), int, double, double, double*, int, double*) const
  5# KEMField::KElectrostaticAnalyticConicSectionIntegrator::ElectricField(KEMField::KConicSection const*, KEMField::KThreeVector_<true> const&) const
****************[KSSTEP WARNING MESSAGE] Failed to execute step <1501> (GSL error: ellint.c:542: domain error (1))

while electrons passing through the hole activate the ksterm_max_z termination. In fact, if I open the root file I find b'gsl_error' and b'term_max_z'. Instead, if I use the anode without a hole like the cathode, I don't get this error, but b'term_death_surface'.

I started to suppose that could be a problem related to the attribute radial_mesh_count of the cylinder_tube_space but it doesn't change and I've already defined the axial_mesh and electrostatic_dirichlet.

The kemfield I'm using is the same as the DipoleTrapSimulation.xml example.

[richeldichel]

Hi @marcogobbo,

thanks for your message. I have also encountered this behavior before, but it was some years ago, so I unfortunately have no steps to reproduce it. In my case, it also terminated the electrons at the right position, so I just treated the "gsl_error"-terminator as the one that I intended to have.
I think though, that it had something to do with the trajectory control. Are you using control_cyclotron or something different? It might make sense to play around with this.

To pinpoint the exact reason of the error, I would like to ask you whether you have a minimal example so that I can reproduce the bug on my machine and be able to properly debug this error.

[marcogobbo]

Thanks for the reply @richeldichel!

I was thinking of using the same technique since the output is what I'm expecting, but I don't really like errors and makes me doubt what I've written. I'm using this type of control

<control_position_error name="control_position_error" absolute_position_error="1e-12" safety_factor="0.75" solver_order="8"/>
<control_length name="stepsizelength" length="1e-4" />
<control_time name="stepsizetime" time="1e-6" />

I posted here the entire structure https://gist.github.com/marcogobbo/e5dd800a79449ac2c0daea3856b2f541

Meanwhile, I'll try to find different approaches using the control tags! Let me know if you find something interesting!

[2xB]

Hey, to add to that: From looking into the source, this error happens when k in KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing is 1 (see https://git.savannah.gnu.org/cgit/gsl.git/tree/specfunc/ellint.c#n539 ).

I'm not quite sure why this happens, but we can investigate by adding code like the following to the beginning of KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing in KEMField/Source/BoundaryIntegrals/Electrostatic/Analytic/src/KElectrostaticAnalyticRingIntegrator.cc directly between double k = sqrt(1. - eta); and double E = E_elliptic(k);:

if (k*k >= 1.0) {
    std::cout << "k*k >= 1.0, info: " << k << ", " << P[0] << ", " << par[0] << ", " << par[1] << ", " << par[2] << ", " << par[3] << ", " << par[4] << ", " << par[5] << ", " << par[6] << std::endl;
}

That should give more information on what exactly goes wrong there.

[marcogobbo]

I've added this code, but it doesn't output anything. For the avoidance of doubt, this is the method we are aiming for

double KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing(const double* P, const double* par)
{
    static KCompleteEllipticIntegral2ndKind E_elliptic;
    static KEllipticEMinusKOverkSquared EK_elliptic;

    double Z = par[2] + P[0] / par[6] * (par[4] - par[2]);
    double R = par[3] + P[0] / par[6] * (par[5] - par[3]);

    double dz = par[0] - Z;
    double dr = par[1] - R;
    double sumr = R + par[1];

    double eta = (dr * dr + dz * dz) / (sumr * sumr + dz * dz);
    double k = sqrt(1. - eta);
    if (k*k >= 1.0) {
    	std::cout << "k*k >= 1.0, info: " << k << ", " << P[0] << ", " << par[0] << ", " << par[1] << ", " << par[2] << ", " << par[3] << ", " << par[4] << ", " << par[5] << ", " << par[6] << std::endl;
	}
    double E = E_elliptic(k);
    double S = sqrt(sumr * sumr + dz * dz);
    double EK = EK_elliptic(k);

    return R / (S * S * S) * (-2. * R * EK + dr / eta * E);
}

The output error is still the same, I've tried to remove the if but still doesn't print anything.

15:35:12 [ERROR] l/Utility/KGslErrorHandler.cxx:84    GSL error: ellint.c:542: domain error (1)
  1# gsl_sf_ellint_Ecomp_e
  2# gsl_sf_ellint_Ecomp
  3# KEMField::KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing(double const*, double const*)
  4# KEMField::KGaussianQuadrature::operator()(double (**)(double const*, double const*), int, double, double, double*, int, double*) const
  5# KEMField::KElectrostaticAnalyticConicSectionIntegrator::ElectricField(KEMField::KConicSection const*, KEMField::KThreeVector_<true> const&) const
****************[KSSTEP WARNING MESSAGE] Failed to execute step <166> (GSL error: ellint.c:542: domain error (1))

[2xB]

Thank you for trying! This is weird. Since that is something that occurred to me multiple times: Are you sure you did not just make (or make -j7 or similar) but also ran make install before trying to run your modified Kassiopeia?

The function not printing anything would be really weird otherwise, since it obviously fails in that function (as from your stack trace: 3# KEMField::KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing)...

[marcogobbo]

I forgot to run one of the make... Forget it

I ran the simulation for 20 events and I got this output

k*k >= 1.0, info: 1, 0.000155797, 0.0085, 0.0148722, 0.0085, 0.015028, 0.0085, 0.0147931, 0.000234907
k*k >= 1.0, info: 1, 0.000126238, 0.0085, 0.00663569, 0.0085, 0.00676193, 0.0085, 0.00649529, 0.000266639
k*k >= 1.0, info: 1, 7.38246e-05, 0.0085, 0.00845964, 0.0085, 0.00853347, 0.0085, 0.00828787, 0.000245593
k*k >= 1.0, info: 1, 0.00021917, 0.0085, 0.0148088, 0.0085, 0.015028, 0.0085, 0.0147931, 0.000234907
k*k >= 1.0, info: 1, 0.000128264, 0.0085, 0.00740961, 0.0085, 0.00753788, 0.0085, 0.00728262, 0.000255259
k*k >= 1.0, info: 1, 2.55747e-05, 0.0085, 0.0179503, 0.0085, 0.0179759, 0.0085, 0.0177174, 0.000258478
k*k >= 1.0, info: 1, 0.000154777, 0.0085, 0.0148732, 0.0085, 0.015028, 0.0085, 0.0147931, 0.000234907
k*k >= 1.0, info: 1, 8.72904e-05, 0.0085, 0.0108175, 0.0085, 0.0109048, 0.0085, 0.0106733, 0.000231569
k*k >= 1.0, info: 1, 8.95896e-05, 0.0085, 0.0195591, 0.0085, 0.0196487, 0.0085, 0.0193445, 0.000304231
k*k >= 1.0, info: 1, 8.53806e-05, 0.0085, 0.0108194, 0.0085, 0.0109048, 0.0085, 0.0106733, 0.000231569
k*k >= 1.0, info: 1, 8.85485e-05, 0.0085, 0.0089306, 0.0085, 0.00901915, 0.0085, 0.00877718, 0.000241973
k*k >= 1.0, info: 1, 0.000129858, 0.0085, 0.0156106, 0.0085, 0.0157405, 0.0085, 0.0155016, 0.000238866
k*k >= 1.0, info: 1, 0.000148901, 0.0085, 0.0163176, 0.0085, 0.0164665, 0.0085, 0.0162228, 0.000243709
k*k >= 1.0, info: 1, 5.07575e-05, 0.0085, 0.0154509, 0.0085, 0.0155016, 0.0085, 0.0152642, 0.000237464
k*k >= 1.0, info: 1, 0.00019155, 0.0085, 0.0132148, 0.0085, 0.0134063, 0.0085, 0.0131785, 0.00022781
k*k >= 1.0, info: 1, 0.000226605, 0.0085, 0.00806127, 0.0085, 0.00828787, 0.0085, 0.00804022, 0.000247651
k*k >= 1.0, info: 1, 0.000164645, 0.0085, 0.0143947, 0.0085, 0.0145594, 0.0085, 0.0143267, 0.000232624
k*k >= 1.0, info: 1, 0.000215249, 0.0085, 0.0194335, 0.0085, 0.0196487, 0.0085, 0.0193445, 0.000304231
k*k >= 1.0, info: 1, 0.000227183, 0.0085, 0.0157537, 0.0085, 0.0159808, 0.0085, 0.0157405, 0.000240364

[2xB]

Thank you! In the meantime I also tested it and an example output where it actually fails later is

k*k >= 1.0, info: 1, 0.000372094, 0.0085, 0.00662791, 0.0085, 0.007, 0.0085, 0.0065, 0.0005
12:31:56 [ERROR] l/Utility/KGslErrorHandler.cxx:84    GSL error: ellint.c:542: domain error (1)
  1# gsl_sf_ellint_Ecomp_e
  2# gsl_sf_ellint_Ecomp
  3# KEMField::KElectrostaticAnalyticRingIntegrator::EFieldRFromChargedRing(double const*, double const*)
  4# KEMField::KGaussianQuadrature::operator()(double (**)(double const*, double const*), int, double, double, double*, int, double*) const
  5# KEMField::KElectrostaticAnalyticConicSectionIntegrator::ElectricField(KEMField::KConicSection const*, KEMField::KThreeVector_<true> const&) const

Anyways, it turns out that this error occurs when one calculates a field e.g. on a charge-defined conic section.
The interesting question now is whether this is intended - in which case we should add a better error message so people can e.g. understand when to use terminators to avoid particles being tracked into walls - or if this is not intended, in which case we of course need to fix this. In case anyone knows the answer, that would be very interesting - maybe @richeldichel ?