Creating non-uniform electric field

[cghw02]

Hello,
I am using Geant4.10.6.p01 on Windows 10.

I am attempting to setup a non-uniform electric field around a sharp tip. The simulation volume is on the order of a few 100 nanometers.

Thus far, I have copied the relevant files from the field02 example and managed to create a uniform electric field in my model. I can fly the electrons and I see the trajectories bend in the field and I can view the field with the /vis/scene/add/electricField command.

Now I try to change the field to a non-uniform electric field. There is no example provided for a non-uniform electric field, but there is one for a non-uniform magnetic field (purging magnet). Hence I have taken this code, but instead of creating a magnetic field, I create an electric field (for the moment using the same file that purging magnet uses). The file appears to be read correctly as I get this output which is very similar to the output of the purging magnet example.

  Electric field

---

—> Reading the field grid from PurgMag3D.TABLE …
[ Number of values x,y,z: 6 6 10 ]

—> … done reading
—> assumed the order: x, y, z, Bx, By, Bz
—> Min values x,y,z: -5e-06 -5e-06 -2.6e-05 cm
—> Max values x,y,z: 5e-06 1.7e-05 1e-05 cm
—> The field will be offset by 0 cm

After reordering if neccesary
—> Min values x,y,z: -5e-06 -5e-06 -2.6e-05 cm
—> Max values x,y,z: 5e-06 1.7e-05 1e-05 cm
—> Dif values x,y,z (range): 1e-05 2.2e-05 3.6e-05 cm in z

F02ElectricFieldSetup: The minimal step is equal to 0.001 mm
G4ClassicalRK4 is called

However, I do not see the electron trajectories curve and when I try to view the electric field I get the message:
No Electric field in this extent.

The electric field constructor (taken from Field02 example) is this (I commented out the line where a uniform field is setup):

F02ElectricFieldSetup::F02ElectricFieldSetup()
: fMinStep(0.0010mm), // minimal step of 1 micrometers
fFieldManager(0),
fChordFinder(0),
fEquation(0),
fEMfield(0),
fElFieldValue(),
fStepper(0),
fIntgrDriver(0),
fStepperType(4), // ClassicalRK4 – the default stepper
fFieldMessenger(nullptr)
{
// fEMfield = new G4UniformElectricField(G4ThreeVector(0.0,-1.0kilovolt/cm,0.0));
G4double zOffset = 0.0 * mm;
fEMfield = new PurgMagTabulatedEField3D(“PurgMag3D.TABLE”, zOffset);
fEquation = new G4EqMagElectricField(fEMfield);

fFieldManager = GetGlobalFieldManager();

UpdateIntegrator();
fFieldMessenger = new F02FieldMessenger(this);
}

The code for PurgMagTabulatedEField3D is essentially the same as PurgMagTabulatedField3D in the purging magnet example.

The major difference is that PurgMagTabulatedEField3D has the following code in the include file.
class PurgMagTabulatedEField3D
#ifndef STANDALONE
: public G4ElectricField
#endif

so PurgMagTabulatedEField3D inherits from G4ElectricField whereas PurgMagTabulatedField3D (from the purging magnet example) inherits from G4MagneticField.

Any advice/suggestions would be welcome.
Thanks
Chris

[allison]

Hi Chris

I take some responsibility for the visualisation of fields. So I am pleased that you can view the field with the /vis/scene/add/electricField command for a uniform electric field.

What worries me is, “No Electric field in this extent.” For the minute I cannot see where that comes from. I cannot see that output in any statement in Geant4.10.6.p01. Ah! In G4VFieldModel::DescribeYourselfTo there is

265   if (FieldMagnitudeMax <= 0.) {
266     G4cout << "No " << fTypeOfField << " field in this extent." << G4end;
267     return;

So it looks like G4VFieldModel thinks there is no field (it’s maximum value is <= 0.). Why might that be?

Are you in a position to debug this? (Compile in Debug mode and use a debugger. Set a breakpoint at line 265, or maybe earlier. Maybe set the breakpoint at the start of G4VFieldModel::DescribeYourselfTo. What does it think the field values are?)

Sorry not to be more helpful regarding the field itself. It might be established correctly and, possibly, the problem is in the visualising modelling code.

But since you do not see the electron trajectories curve, maybe the problem lies elsewhere.

John

[cghw02]

Hi John,
Thanks for getting back so soon. I am trying to debug using Windbg. I am not so familiar with this debugger and have mostly used gdb on Linux before. Anyway, immediately prior to the statement:
G4cout << “No " << fTypeOfField << " field in this extent.” << G4endl;

it is going round a loop determining the field at certain points. It would seem the magnitude of the field is zero at all points. I am wondering whether the field is defined on a much larger scale (cm) as compared to my model (nm). However, upon checking my class PurgMagTabulatedEField3D.cc, the distances over which the field is defined is given by:

maxx = xval * lenUnit;
maxy = yval * lenUnit;
maxz = zval * lenUnit;

where lenUnit in my case is 0.001mm whereas for the original class PurgMagTabulatedField3D.cc lenUnit was defined as 10cm.

I should say that I am taking this afternoon and tomorrow as holiday, so there will be a hiaitus in replying to any new posts for a little while.
Best Wishes
Chris

[cghw02]

Hi John,
I forgot to mention that I was able to check certain values in G4VFieldModel, which are:

xHalfScene 17.226397160708160072
yHalfScene 24.050794322416322757
zHalfscene -6.5757943204163211703
xSceneCenter -7.2263971607081600723
ySceneCenter 4.0507943204163225914
zSceneCenter 2.0253971602081612957

Are these numbers in meters or just relative to some scene size ? For other numbers, it often says value unavailable or just provides the address, not the number contained in the address…

[allison]

Hi Chris

These numbers xHalfScene, etc., will be in Geant4 internal units, which (if you have not changed them), are mm. (Whenever you print something you should divide by the unit, e.g:
G4cout << "xHalfScene/mm " << xHalfScene/mm << G4endl;

So…it looks like G4VFieldModel thinks the extent is a few tens of mm. But you say your field is only a few nm. I think G4VFieldModel assumes the extent is that of the solid in the logical volume in which the field is defined.

Let’s pick this up when you return from your break.

John

[cghw02]

Hi John,
I have tried to change the size of my model - just to try and get past the “No Electric field in this extent.” problem. What I had been finding was that (as you say) the points in my model lay outside the minimum or maximum extent in some dimension, so all electric field points were returned as zero, which then leads to the “No Electric field in this extent.” message.

However, now I have put these lines in
void PurgMagTabulatedEField3D::GetFieldValue(const double point[4],
double *Efield ) const

#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "minx: " << minx << " maxx: " << maxx << " x: " << x << " Efield[0]: " << Efield[0] << G4endl;
G4cout << "miny: " << miny << " maxy: " << maxy << " y: " << y << " Efield[1]: " << Efield[1] << G4endl;
G4cout << "minz: " << minz << " maxz: " << maxz << " z: " << z << " Efield[2]: " << Efield[2] << G4endl;
#endif

An example of the last few lines of the output is below after issuing the command:
/vis/scene/add/electricField

minx: -5000 maxx: 5000 x: 1000 Efield[0]: 0.236851
miny: -5000 maxy: 17000 y: 435.954 Efield[1]: -0.00305162
minz: -26000 maxz: 10000 z: 1398.45 Efield[2]: -0.0398415
No Electric field in this extent.
Idle>

and so a non zero electric field would appear to be returned at least some of the time - but I still get the same message, so I am still rather unsure what is happening.

Chris

[cghw02]

Hi John,
I have discovered this bit of code:
void G4ElectricFieldModel::
GetFieldAtLocation(const G4Field* field, const G4Point3D& position,
G4double time, G4Point3D& result) const {
if (!field) return; // No action if no field

G4double xyzt[4] = { position.x(), position.y(), position.z(), time };
G4double BEvals[6] = {0.}; // Field returns {Bx,By,Bz,Ex,Ey,Ez}
field->GetFieldValue(xyzt, BEvals);

result.set(BEvals[3], BEvals[4], BEvals[5]);
return;

So it seems that in my GetFieldValue routine, I need to put the electric field values in the locations 3, 4 and 5 - not 0, 1 and 2 in a 6 element vector. I will try that.

Best Wishes
Chris

[cghw02]

Dear John,
I am now seeing the electric field and I have modified the code so that it now reads a file which represents the electric field around a sharp tip near a conducting surface. The field was generated using COMSOL. I have tried to make sure that the electric field and the geant4 model are centred around the coordinates 0,0,0, but so far the field displayed in Geant4 is not like that shown in COMSOL. The numbers appear to be read in correctly. Is there a recommended approach to checking the coordinates that the model and the electric field are using ?
Thanks
Chris

[allison]

HI Chris
Good other you have found the field!
Visualisation is there to help, but for complete and precise checking all I can recommend is the old-fashioned time-tested method of adding print statements or using a debugger.
But I’m curious about “the field displayed in Geant4 is not like that shown in COMSOL”. In what way is it different?
John

[allison]

That should read “Good that you have found the field”.

[mkelsey]

COMSOL makes beautiful pictures, typically showing things like field magnitude, or the potential, in the form of a continuous rainbow filling the volume of interest.

The G4 field visualization are arrows of varying lengths, drawn in a very sparse grid spanning the entire world volume.

Interpreting the two in a consistent way is an exercise in difficulty and frustruation

[allison]

Hi Mike

It was my suspicion that there is a significant difference in sophistication between G4 and COMSOL, and I can appreciate that interpreting the two in a consistent way is an exercise in difficulty and frustration. I’ve never tried it myself, but it should be possible to see a correspondence:

1. Note that the arrow represents the field magnitude and direction at its base (not at the head);

2. The grid on which it is drawn does not have to be sparse. See guidance on /vis/scene/add/electricField.

3. You don’t have to draw the field in the whole world volume. We now have:

/vis/set/extentForField
/vis/set/volumeForField
/vis/touchable/extentForField
/vis/touchable/volumeForField
(For the last two you have to /vis/set/touchable.)

It would be good - Chris and Mike - to hear if those new features are useful.

Cheers
John

[allison]

Hi Chris

After all’s said and done on the vis side, you need to check that tracking is seeing your field. I’m not an expert here but I quote from an email of a colleague who has been following this thread.

"I saw the item in the forum on “non-uniform electric field” and
your answers. I know little about how fields are navigated but
from experience with G4Beamline I can speculate that if his
field extent is really only a few nanometers he is going to have
problems in the tracking. The stepper will not see the field and
very likely will just step over it. I have seen this even with
fields of a few mm. The stepper may be very clever about sampling
field values, but this is like finding a needle in a haystack.

“In G4Beamline this can be usually be dealt with e.g. by putting
a pseudo volume around the field region thus forcing the tracking
to end the step at each boundary, or by just setting a very small
step limit. This works on the millimeter scale but I don’t know
about nanometers! One thing I haven’t tried is to make the
field extent much larger and fill the extra with very small nonzero
values.”

[cghw02]

Hi John and Mike,
Thanks for your comments (and those of your colleague John). I have uploaded a couple of pictures from COMSOL and geant4 simulations as they currently stand. I have written a bit of code for the field messenger so when I type (for example)
/field/getField 10 10 10 nm,
it will return the field at the given location and the field would appear to be the same as that in COMSOL. However, although the field is currently there in Geant4, only a couple of arrows are shown. I guess there is a scaling problem, but so far not tracked it down.
Thanks again
Chris

[allison]

Hi Mike
Good to know that your /field/getField returns what you expect.
Can you not increase the sampling with /vis/scene/add/electricField 100 lightArrow (lightArrow for speed, perhaps)? Or /vis/scene/add/electricField 1000 or /vis/scene/add/electricField 10000, etc.
Let me know if you feel we may need to offer more options.
J

[cghw02]

Hi John,
I tried what you suggested and I see many more arrows and they appear clustered around the tip.
I´m not sure what the different colours mean, but it looks good. However, now I find that the electrons are not curving in this field. In addition, I try to launch the electrons from a particular position, but I do not see the track from this position and only after there is some interaction in the sample. I think maybe as the material properties were defined in the GDML file previously, then some bug has crept in as I now define the material properties within the program. I have also tried increasing the size of the simulation (up to mm), but this does not seem to change these problems. I will continue to check what might be wrong.
Best Wishes
Chris

[allison]

Hi Chris

Nice pictures! Do the field arrows make sense? (They all seem to be pointing one way - that can’t be right. You can’t have a non-uniform uni-directional field. It has to have some divergence - Maxwell’s equations.)

Also, if it gets cluttered you can use lines (option lightArrow) - see guidance (help /vis/scene/add/electricField or the Help tab on the Qt GUI). Also the guidance explains the colours - except it’s wrong (mea culpa) - the colours map field strength as a fraction of the maximum magnitude: 0->0.5->1 is red->green->blue (not blue->green->red, as it says in the guidance - I’ll fix it). But in principle, it really pays to read the command guidance.

As for tracking/transportation, I hand over to the experts.

Good luck
John

[cghw02]

Dear John,
I hope the attached image makes more sense in terms of Maxwell´s equations ! I think perhaps the view was from too far away before.
I currently find that if I set the gun position/direction/energy/particle and then do
run/beamOn 1
I might get a trajectory shown, but if I do again run/beamOn 1 , the program crashes. If I do a backtrace while in gdb (running on Linux), it seems to crash when calling the destructor for G4SmoothTrajectoryPoint. This would appear to be part of the visualisation code. Any idea what might be happening ? It seems to occur whether I am displaying the electric field or not.
Thanks
Chris

[cghw02]

Dear John,
Further to my last message, I am now seeing electron trajectories bend in my non-uniform field and bounce on the surface (as they should). With regard to the crashing, I still see this, but it seems to happen after I increase the energy from the inital 10 eV to 100 eV.

Best Wishes
Chris

[allison]

Hi Chris

Nice picture, but, as Mike said, not easy to see if the field is sensible.

I’d need a full traceback. G4SmoothTrajectoryPoint is not part of vis but is used by vis. I have no idea what might be happening. Are there any error/warning messages. I can use /vis/scene/add/trajectory smooth without crashing, but, of course, I have not tried in a non-uniform electric field.

John