Message: Re: Problem with magnetic field from a table Not Logged In (login)
 Next-in-Thread Next-in-Thread
 Next-in-Forum Next-in-Forum

None Re: Problem with magnetic field from a table 

Forum: Fields: Magnetic and Otherwise
Re: None Problem with magnetic field from a table (Mariusz)
Date: 25 Jan, 2008
From: John Apostolakis <John Apostolakis>

Dear Mariusz,

Geant4 uses the GetFieldValue method that you create many times during 
every step.  Looking at the method that you have written, I can see that 
it is making a very large number of operations, in order to calculate 
the value.  So it is not surprising that it your program is very slow.

I have two suggestion to speed up your key GetFieldValue method:

1) Avoid large loops:

        for(G4int i=0; i<4900; i++) {
                if(fabs(cPoint[0]-x[i])<dmesh/2 ...

     Use the functions for converting a floating point number to a 
integer, in order to calculate the nearest 'i' integer to your point.  
Some code along the lines of:
                 double ratio = (cPoint[0]/dmesh);
                 int i = round( ratio );   // Double check this - it 
should round to nearest integer maybe it is std::round
                 if( fabs( cpoint[0]- x[i]) < dmesh/2 ) ...
etc

2)  to replace the line

G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);

with

G4double radius = std::sqrt(Point[0]*cPoint[0]+cPoint[1]*cPoint[1]);

since
  - the power function is very costly compared to a multiplication, and
  - square root is faster than pow( .. , 0.5)
The power method generally requires a logarithm and an exponential.

I expect that this will speed up your field method, and your 
application, by a factor of 100 or more.  I expect that your particles 
will actually move.

Best regards,
John

Mariusz wrote:
> *** Discussion title: Fields: Magnetic and Otherwise
> Email replies to PublicHyperNews@slac.stanford.edu must include:
>   In-Reply-To: <"/emfields/146"@geant4-hn.slac.stanford.edu>
>   Subject: ...change this to be about your reply.
>
>  Hi,
>
>   I have a 2d table with values of Bx and By (Bz is 0). It is quite large: 4900 points
> with spacing of 7.7mm. I use Geant 4.9.0.p01.
>   I have a class:
>
> ---------------------VVV ------------------------
> class WSMagneticField: public G4MagneticField
> {
>   public:
>
>    WSMagneticField(G4ThreeVector);  //  The value of the field
>    WSMagneticField();               //  A zero field
>   ~WSMagneticField();  
>
>    void SetFieldValue(G4double fieldValue);
>    void SetFieldValue(G4ThreeVector fieldVector);
>
>    void GetFieldValue(const double Point[3],
>                 double *Bfield) const;           
>
>   protected:
>       // Find the global Field Manager
>       G4FieldManager* GetGlobalFieldManager();   
>   private:
>       static const int imesh=4900;
>       float x[imesh];
>       float y[imesh];
>       float Bx[imesh];
>       float By[imesh];   
> };
> ----------------------AAA---------------------
> In the constructor I fill the tables x,y,Bx,By:
>
> ---------------------VVV----------------------
> WSMagneticField::WSMagneticField() 
> {
>         // here open and read file      
>         ifstream bmap;
>         string line;
>         const char *charline;
>         int ni1,ni2,idx=0;
>         float modB;
>         char ni3;
>
>         bmap.open("Bfield.data");
>         while(bmap) {
>                 getline(bmap,line);
>                 charline=line.c_str();
>                 sscanf(charline,"  %d  %d  %f  %f  %f  %f  %f  %c\n",
>                    &ni1,&ni2,&x[idx],&y[idx],&Bx[idx],&By[idx],&modB,&ni3);
>                 G4cout << "mag field check i = " << idx << " x = " << x[idx] << " Bx = " << Bx[idx] <<G4endl;
>                 idx++;                   
>         }       
>
>   GetGlobalFieldManager()->SetDetectorField(this);
>   GetGlobalFieldManager()->CreateChordFinder(this);
> }
> --------------------AAA-------------------------
>
> and I have GetFieldValue function:
> --------------------VVV-------------------------
>
> // argument is position in mm
> void WSMagneticField::GetFieldValue(const double Point[3], double *Bfield) const
> {
>         Bfield[0] = 0.;
>         Bfield[1] = 0.;
>         Bfield[2] = 0.;
>         G4float dmesh=7.93; // [mm]
>         G4double cPoint[3];
>         cPoint[0]=Point[0];
>         cPoint[1]=Point[1]-77.; // [mm] shift back the Quad to World ref.    
>         cPoint[2]=Point[2];     
>         G4double radius = pow(pow(cPoint[0],2)+pow(cPoint[1],2),0.5);
>         if(cPoint[0]>0 && radius < 240. && fabs(cPoint[2]) < 1400.) {
>          for(G4int i=0; i<3; i++) if(cPoint[i]<0) cPoint[i]*=(-1); // because of symmetry
>          for(G4int i=0; i<4900; i++) {
>
>                 if(fabs(cPoint[0]-x[i])<dmesh/2 && fabs(cPoint[1]-y[i])<dmesh/2) {
>                                 Bfield[0]=Bx[i];
>                                 Bfield[1]=By[i];
>                                 break;
>                 }
>          }
>         }       
> }
> -------------------AAA------------------------- 
>
> And then, in DetectorConstruction I have:
> -------------------VVV-------------------------
>         WSMagneticField* myField = new WSMagneticField();
>         G4FieldManager* fieldMgr
>                                 = G4TransportationManager::GetTransportationManager()
>                                 ->GetFieldManager();
>
> // I tried this and it did not work either //
> fieldMgr->SetDetectorField(myField); //
> fieldMgr->CreateChordFinder(myField);
>
> /* ... geometry declaration ... */
>
>    G4FieldManager *localFieldMgr = new G4FieldManager(myField);
>    logicQ5->SetFieldManager(localFieldMgr,true);
> --------------------AAA-------------------------
>
> The program compiles, the file with table of values is read but never
> ever any event has finished being calculated. I've left it even for the
> whole night even if in normal conditions one event takes about 15
> minutes.
>
> I would appreciate any suggestion leading to solution;-)
>
>  Mariusz Sapinski  
>
>   

Inline Depth:
 1 1
 All All
Outline Depth:
 1 1
 2 2
 All All
Add message: (add)

1 None: Re: Problem with magnetic field from a table   (Mariusz - 26 Jan, 2008)
1 None: Re: Problem with magnetic field from a table   (John Apostolakis - 28 Jan, 2008)
3 None: Re: Problem with magnetic field from a table   (Mariusz - 27 Jan, 2008)
1 None: Re: Problem with magnetic field from a table   (John Apostolakis - 28 Jan, 2008)
 Add Message Add Message
to: "Re: Problem with magnetic field from a table"

 Subscribe Subscribe

This site runs SLAC HyperNews version 1.11-slac-98, derived from the original HyperNews


[ Geant 4 Home | Geant 4 HyperNews | Search | Request New Forum | Feedback ]