The OPAL Discussion Forum

[masao nakao <nakaom AT rcnp.osaka-u.ac.jp>]

Dear Achim,

Thank you for your kindly reply about SAAMG and vtk2h5grid.
I have installed vtk and vtk2h5grid and converted mesh file of inflector 
to h5 file.
Then, I found that OPAL returns error such as:

Error> *** User error detected by function 
"ArbitraryDomain::ArbitraryDomain()"
Error>     This domain is currently not available.
Error>     This domain is currently not available.
application called MPI_Abort(MPI_COMM_WORLD, -100) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=-100
:
system msg for write_line failure : Bad file descriptor

I have attached the OPAL input file to this email.

Yours  sincerely,
Masao






On 2020/08/18 17:15, Gsell Achim (PSI) wrote:
> Dear Masao
>
> > On 18 Aug 2020, at 08:10, masao nakao <nakaom AT rcnp.osaka-u.ac.jp 
> > <mailto:nakaom AT rcnp.osaka-u.ac.jp>> wrote:
> >
> > Dear Achim,
> >
> > I would like to simulate an AVF cyclotron injection system using 
> > inflector electrodes. And so far I have been able to simulate the 
> > case where a single particle can be injected into the central plane 
> > with the electric field calculated by TOSCA.
> >
> > Next, I would like to calculate including the space-charge effect at 
> > an injection current of ~10 mA. I was able to save the mesh shape of 
> > the inflector electrodes as a vtk file, but I have to convert it to 
> > h5block format and input it to OPAL. Can you tell me how to convert 
> > vtk to h5block format and input it to OPAL?
>
> H5hut includes a converter. Unfortunately you have to compile it 
> yourself. Here is a short recipe to compile H5hut with the converter:
>
>   * compile/install a serial(!) version of HDF5 1.8 or 1.10
>     (https://www.hdfgroup.org/downloads/hdf5/)
>   * compile/install VTK 8.2 (https://vtk.org/download/)
>   * clone the H5hut repository (https://gitlab.psi.ch/H5hut/src.git)
>   * configure H5hut with (replace *_PREFIX with the installation
>     directories of HDF5, VTK)
>
> configure --enable-vtkconverter --with-hdf5=HDF5_PREFIX 
> --with-vtk=VTK_PREFIX
>
>   * compile & install
>
> make && make install
>
>   * convert your (legacy!) VTK file
>
> vtk2h5grid FILENAME
>
> If you have problems compiling the converter, I can convert the file 
> for you.
>
> Best
>
> Achim
>
> > Thank you.
> >
> >
> > Yours, sincerely,
> > Masao
> >
> >
> >
> >
> > On 2020/03/20 21:41, Daniel Winklehner wrote:
> > > Dear Masao,
> > >
> > > Yes, you can do it, but it requires some fiddling around at the moment.
> > >
> > > The most important question is: what beam current are you 
> > > simulating? Is space charge an important factor or not so much?
> > >
> > > If space charge is very important for your application, you will 
> > > have to compile OPAL with the SAAMG option ON. To my knowledge, 
> > > there is no binary package that supports usage of the SAAMG solver 
> > > that you can just download, but Achim can correct me on that.
> > >
> > > If space-charge is unimportant, you can use the OPAL 2.2 binary.
> > >
> > > In both cases, you will need to prepare the electric AND magnetic 
> > > fields as 3D fieldmaps in the center of the cyclotron as *.h5part 
> > > format files. OPAL-cycl, by default, uses the midplane field only 
> > > and gets the z =/= 0 components from a series expansion. For the 
> > > electrostatic inflector, you will need a more accurate field 
> > > description along the spiral trajectory. There have been some email 
> > > list discussions recently on how to transform 3D field data from 
> > > ANSYS/TOSCA/COMSOL to h5part using a C script called 
> > > ascii2h5part.cpp. I am using python with h5py, which works too. See 
> > > attached for electrostatic case. Resulting fields can be checked 
> > > with VisIt 3.1.1.
> > >
> > > The fields can be loaded into OPAL-cycl by using the BANDRF type 
> > > cyclotron and loading the fields as RF cavity fields with frequency 
> > > = 1e-36 (i.e. static fields) and RFPHI = 0.0 and PI/2 for electric 
> > > and magnetic fields, respectively. The units for these fields should 
> > > be kV/mm (MV/m) and mm.
> > >
> > > Finally, the geometry can be loaded as well. Firstly, as the 
> > > boundary for particle termination, if you are interested in the 
> > > losses. Secondly, as the boundary conditions for the SAAMG field 
> > > solver (if you are using it). For geometry preparation, you will 
> > > need the "negative" or "vacuum" space rather than the individual 
> > > electrodes. This is described in the paper you mentioned. Just 
> > > create a cylinder that encloses all your electrodes (maybe start 
> > > with just the spirals) and subtract the electrodes from the cylinder 
> > > as a Boolean operation in the program of your choice. Inventor 
> > > works, so does FreeCAD (which is open source). Save as iges, stp, 
> > > etc. and open in GMSH. In gmsh, you can mesh your geometry and save 
> > > the mesh as vtk. Achim write a C++ program that translates vtk to 
> > > h5block, which is the format OPAL needs for geometry.
> > >
> > > If you do use SAAMG, you have to use the "SPIRAL" cyclotron type 
> > > instead of BANDRF.
> > >
> > > I will write up better instructions for the whole process over the 
> > > weekend and provide a simple example input deck to help along with 
> > > getting started on this non-trivial process.
> > >
> > > Best,
> > >    Daniel
> > >
> > >
> > >
> > > -----Original Message-----
> > > From: opal-request AT lists.psi.ch <mailto:opal-request AT lists.psi.ch> 
> > > <opal-request AT lists.psi.ch <mailto:opal-request AT lists.psi.ch>> On 
> > > Behalf Of masao nakao
> > > Sent: Friday, March 20, 2020 07:44
> > > To: opal AT lists.psi.ch <mailto:opal AT lists.psi.ch>
> > > Cc: nakaom AT rcnp.osaka-u.ac.jp <mailto:nakaom AT rcnp.osaka-u.ac.jp>
> > > Subject: [Opal] Fwd: 3DElectroStatic
> > >
> > > Dear OPAL developers,
> > >
> > > I would like to simulate injection system of AVF cycltoron with 
> > > electrostatic inflector.
> > > Though OPAL manual shows 3DElectroStatic is "Not implemented 
> > > yet."(p.173), a paper [1] shows results by OPAL. How can I calcurate 
> > > like that? I have 3D electric field data culculated by TOSCA.
> > > Thank you.
> > >
> > > [1] "Realistic simulations of a cyclotron spiral inflector within a 
> > > particle-in-cell framework"
> > > Daniel Winklehner et al.
> > > PHYSICAL REVIEW ACCELERATORS AND BEAMS 20, 124201 (2017)
> > >
> > > Yours, sincerely,
> > > Masao
> > >
> > >
> > > --
> > > Masao Nakao
> > > Osaka University
> > > Research Center for Nuclear Physics
> > > Mail : nakaom AT rcnp.osaka-u.ac.jp <mailto:nakaom AT rcnp.osaka-u.ac.jp>
> > > Phone: 06-6879-8914
> > > FAX  : 06-6879-8899
>
> Senior Scientist
> Paul Scherrer Institut
> Laboratory for Scientific Computing and Modeling
> OHSA/D04
> 5232 Villigen PSI


--
中尾 政夫 (Nakao Masao)
群馬大学 重粒子線医学研究センター
Mail : nakaom AT gunma-u.ac.jp
電話 : 027-220-8378
FAX  : 027-220-8379
携帯 : 070-7774-6423
内線 : 23755

// 2d magneticfield + 3d inflector + 3d dee field + geometry

OPTION, ECHO= FALSE;
OPTION, PSDUMPFREQ= 720;
OPTION, PSDUMPEACHTURN= TRUE;
OPTION, STATDUMPFREQ = 1;
OPTION, VERSION = 20000;
TITLE,STRING= "AVF";

real Edes = 0.000050; //GeV
real r0   = 0.049888;         //m
real phi0 = 0.0;    // initial azimuth  [-180, 180)

real f1=  16.845352; // experimental value

real gamma=  (Edes+PMASS)/PMASS;
real beta=   sqrt(1-(1/gamma^2)); 
real gambet= gamma*beta; 
real P0 =    gamma*beta*PMASS;
real brho = P0 / 0.29979246;
value,{Edes,brho,P0};

// 90deg dee
real r1= 35.0;
real r2= 1050.0;
//real f2=0.000000000001;
real f2=f1 * 2.0;
real volt01= 0.050;
real volt02= 0.050;

real phiv = (PI / 180.0) * 95 ;


Ring: CYCLOTRON, TYPE="BANDRF", CYHARMON=2, 
        PHIINIT=0, ZINIT=0.0, RINIT=0.0 ,  RFFREQ = 1e-36, 
        PRINIT=0, PZINIT = -gambet,
        FMAPFN="p65-580_0mm_opal.table", SYMMETRY=1.0,
        RFPHI={0.0,phiv}, 
RFMAPFN={"inf90b498284020v4-z59d9-rot25.table3.h5part","dee90-243in-0-opal.tabledee.h5part"},
  SPIRAL=true,
        ESCALE={7.505E-3,5.0E-2},  // in MV
        SUPERPOSE={true,true};
        

// Collimators
COL1a0: CCollimator, XSTART=0.0, XEND=0.0,YSTART=1100.0, 
YEND=-1100,ZSTART=25.0, ZEND=60.0, WIDTH=2000.0;
COL2a0: CCollimator, XSTART=0.0, XEND=0.0,YSTART=1100.0, 
YEND=-1100,ZSTART=-60.0, ZEND=-25.0, WIDTH=2000.0;
COL3a0: CCollimator, XSTART=0.0, XEND=0.0,YSTART=-1000.0, 
YEND=-1100,ZSTART=-50.0, ZEND=50.0, WIDTH=10.0;

PR22ax0: Probe,  XSTART=0.0, XEND=1100.0,YSTART= 0.0, YEND= 0.0;
PR22ay0: Probe,  XSTART=0.0, XEND=0.0,YSTART= 0.0, YEND= -1100.0;


L1:   LINE = (Ring,COL1a0,COL2a0,COL3a0,PR22ax0,PR22ay0);

// L1:   LINE = (Ring);

// gauss 
Dist1:DISTRIBUTION, TYPE=Gauss, INPUTMOUNITS = NONE,
sigmax=3.0E-3, sigmapx= gambet * 1.0E-2,
sigmay=3.0E-3, sigmapy= gambet * 1.0E-2,
sigmaz=3.0E-3, sigmapz= gambet * 1.0E-2,
CORRX = 0,  CORRY =0,  CORRT = 0,
OFFSETZ=70.0*1.0E-3;


Fs1:FIELDSOLVER, FSTYPE=SAAMG, MX=16, MY=16, MT=16, PARFFTX=true, 
PARFFTY=true, PARFFTT=true,
            BCFFTX=open, BCFFTY=open, BCFFTT=open, GEOMETRY=Geo1;

Geo1:GEOMETRY, FGEOM="inf-void-rcnpavf3.h5", INSIDEPOINT = {0.001,0.001,30.3};

Beam1: BEAM, PARTICLE=PROTON, PC=p0, NPART=5000, BCURRENT=1.0E-9, CHARGE=1.0, 
BFREQ= f1*2.0;

SELECT, LINE=L1;

TRACK, LINE=L1, BEAM= Beam1, MAXSTEPS= 5*720, STEPSPERTURN= 720, 
TIMEINTEGRATOR="RK-4";
// RUN, METHOD= "CYCLOTRON-T", BEAM= Beam1, FIELDSOLVER= Fs1, DISTRIBUTION= 
Dist1;
 RUN, METHOD= "CYCLOTRON-T", BEAM= Beam1, FIELDSOLVER= Fs1, DISTRIBUTION= 
Dist1,TURNS=1, BOUNDARYGEOMETRY = Geo1;
ENDTRACK;
STOP;