RE: [fluka-discuss]: Importing external magnetic field map in MGFLD

From: Niculescu, Gabriel - niculegx <niculegx_at_jmu.edu>
Date: Tue, 11 Aug 2020 21:57:42 +0000

Dear Tamer-
While I am not a FLUKA expert by any means I am familiar enough with it to try to shed some light on your problem. With apologies if the answer sounds pedantic and/or is too long:

* the purpose of the magfld subroutine is two fold (there is a third possible way if you know the field as an analytic function which is not the case here):

part 1: to read magnetic field information produced by other (exterior to FLUKA) programs.

This is the part governed by the "if (FIRST) then " statement.

part 2: given a point (X,Y,Z) the subroutine is supposed to return the magnitude of the field and its directory cosines: BTX, BTY, BTZ, B (arguments 4-7 of magfld).
Note that Fortran is pass-by-reference so the value of the arguments can/will change.

Other than these two requirements (and the need for the directory cosines to sum. in quadrature, to 1 - meaning you only have two independent ones) what you implement in magfld is up to you.

If the 3D arrays for BX, BY, BZ do not suit your problem you can easily make them unidimensional and make the X,y,z coordinates the same size (20059?). I think you can change the size/shape of these to suit your needs.

That said you need to answer, and implement, the logic that will identify which mesh point (points if you want interpolation) that is closest to your current spatial point.


* If you follow the FLUKA prescription of having a uniform mesh then you can easily identify the correct index. For one dimension if, say, your magnetic region extends from
x1 to x2 with a step size xstep then your current index, for this axis, would be:
ix = (current_x - x1)/xstep
- you might want to round that up/down depending on which side of the interval you like best.
- you might also want to wrap that statement in an if block that ensures current_x is in (x1,x2) to start with, thus eliminating the need to compute the field where it should not exist
- repeat the procedure for the other two coordinates to find their indices. the BX(ix,iy,iz) is your magnetic field component on x, ditto for BY and BZ and you can proceed to calculate B and the directory cosines
- see code snippet at the end of this message

* If, on the other hand, you have an irregular mesh (even sorted according to some criteria) you will have to search through that list to find out which node is closest to your current point. Not knowing anything about the problem you are trying to solve I would say this will require figuring out the distance between your current point and (all) points in the mesh. It seems to me that this version will require substantially more operations than the former. One can argue that it is a classic example of trading memory (size of the map) for speed (how long you search though said map).

Remember that you will have to do the above (whichever version or variation thereof) many-many times depending how many particles you need to simulate, step sizes, interactions, etc. So spending time to think/implement a faster solution will pay huge dividends down the road. If your commercial software cannot produce an uniform grid perhaps you can go through an intermediate step of writing a short python/mathematica/matlab/C++/whatever program that can take the comsol mesh and produce an (interpolated) uniform grid.

Again, apologies for the length of the message. Below you will find a copy of magfld that I've used in some of my models (I wish I could take credit for it but that is not the case. someone - forgot name, sorry! gave me the original. I just modified it to suit purpose) ; please bear in mind that both the physical dimensions of the field region, the dimensions of the mesh, and the order in which the mesh is read are "hard-coded" in the file.

Sincerely,
Gabriel Niculescu
Professor of Physics
Department of Physics and Astronomy
James Madison University


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*$ CREATE MAGFLD.FOR
*COPY MAGFLD
*
*===magfld=============================================================*
*
      SUBROUTINE MAGFLD ( X, Y, Z, BTX, BTY, BTZ, B, NREG, IDISC )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'

      LOGICAL LFIRST
      DOUBLE PRECISION DUMMY1, DUMMY2, DUMMY3
      CHARACTER BUFF
      INTEGER XI , YI , ZI , I, J, K, U
      PARAMETER (NX =61)
      PARAMETER (NY =41)
      PARAMETER (NZ =61)
      PARAMETER (U =98)
      DIMENSION BX(NX ,NY ,NZ)
      DIMENSION BY(NX ,NY ,NZ)
      DIMENSION BZ(NX ,NY ,NZ)
      SAVE LFIRST
      DATA LFIRST /. TRUE ./
      IDISC = 0
      IF ( LFIRST ) THEN
         Write(*,*)'In magfld, before reading field map'
      CALL OAUXFI ('mag.txt',U,'OLD',IERR )
      DO 30 I=1 ,61
      DO 20 J=1 ,41
      DO 10 K=1 ,61
* read and discard x,y,z coordinates as we know the size of mesh
* read bx, by, bz values
      READ (U ,*) DUMMY1 , DUMMY2 , DUMMY3 , BX(I,J,K),
     & BY(I,J,K), BZ(I,J,K)
*gauss to tesla conversion
      BX(I,J,K) = BX(I,J,K )*1.0E -4

      BY(I,J,K) = BY(I,J,K )*1.0E -4
      BZ(I,J,K) = BZ(I,J,K )*1.0E -4
10 CONTINUE
20 CONTINUE
30 CONTINUE
      CLOSE (U)
* keep user informed
       Write(*,*)'In magfld, after reading field map'
       LFIRST = . FALSE .
      ENDIF
*
* x in (3,3)
* y in (-5,5)
* z in (-30,30)
* code below needs to match the order in
* which the mesh was read!!! double check this!!!
      XI = NINT (X*10+31)
      YI = NINT (Y*4+21)
      ZI = NINT (Z+31)
    
           IF (X .LT. 3 .AND . X .GT. -3 .AND. Y .LT. 5 .AND.
     & Y .GT. -5 .AND. Z .LT. 30 .AND. Z .GT. -30) THEN
       B = DBLE ( SQRT (BX(XI ,YI ,ZI )**2 + BY(XI ,YI ,ZI )**2
     & + BZ(XI ,YI ,ZI )**2))
*Normalize
      BTX = BX(XI ,YI ,ZI )/B
      BTY = BY(XI ,YI ,ZI )/B
      BTZ = BZ(XI ,YI ,ZI )/B
      ELSE
      B = ZERZER
      BTX = ZERZER
      BTY = ZERZER
      BTZ = ZERZER
      ENDIF
      RETURN




~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
________________________________________
From: owner-fluka-discuss_at_mi.infn.it [owner-fluka-discuss_at_mi.infn.it] on behalf of Tamer Tolba [tamer.tolba_at_uni-hamburg.de]
Sent: Tuesday, August 11, 2020 2:49 PM
To: fluka-discuss
Subject: [fluka-discuss]: Importing external magnetic field map in MGFLD

Dear FLUKA experts,

I am trying to import an external magnetic field map created by FEA software (COMSOL), with data points sorted in mesh rather than coordinate style* and saved in txt file in fluka MGFLD routine to apply it in my FLUKA calculations.

(* i.e. the mag. filed components are sorted in the text file depending on the exported mesh node coordinate, i.e. not following a specific coordinate order, e.g. (X1,Y1,Z1), (X2, Y1, Z1)....).

I have 20059 data points. The data are sorted in the txt file as X(cm) Y(cm) Z(cm) Bx(T) By(T) Bz(T). Since I am not a fortran expert, I tried to follow the recommendation from an old fluka case here<https://urldefense.proofpoint.com/v2/url?u=http-3A__www.fluka.org_web-5Farchive_earchive_new-2Dfluka-2Ddiscuss_2419.html&d=DwMDaQ&c=eLbWYnpnzycBCgmb7vCI4uqNEB9RSjOdn_5nBEmmeq0&r=28hhbdidO9kO7x107as0Eg&m=PAA7e7HspZt9P5N9sSY1o1-725RuOeFDaDXbhoqfOQ4&s=EBRsPYPwFXb-1fjpIp1QjaFy9wZg6SUN1flTpAb8uns&e=>. Attached, please find my modified (from the fluka old case example) magfld_n file. I tag with "TT" wherever I made an entry.

My questions are:

1- From the old fluka example, it seems that the code reads only the values of the mag. filed components, rather than the coordinate and the filed components values, am I right? If yes, does that mean that fluka assumes "by default" a specific coordinate order, which the B_component should be following during exportation by the COMSOL program?

2- If question #1 is correct, is there any way to make fluka reads the data file as it is sorted (whether on fluka or fortran levels)? Which means that it will need to read the all six variables, 3 coordinate values and 3 B_component value.

3- Can a fluka expert have a glance on the code and tell me if I have any mistakes. As I said I am not a fortran expert, so I might also have mistakes on the programing level. I would really appreciate if you can guide me to any exiting example code for my problem, which seems was an old topic. At least this one from 2009 must have solved it somehow ;).

With my best regards,

Tamer Tolba

--
Dr. Tamer Tolba
Institut für Experimentalphysik,
Universität Hamburg ,
Luruper Chaussee 149,
22761 Hamburg - Germany
Tel.: +49 (0)40-8998-4872
__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info
Received on Wed Aug 12 2020 - 01:43:51 CEST

This archive was generated by hypermail 2.3.0 : Wed Aug 12 2020 - 01:44:05 CEST