> Date: Thu, 27 Jun 2013 09:48:04 +0200
> To: gadget-list_at_MPA-Garching.MPG.DE
> Subject: Re: [gadget-list] IC files_blocks
>
>
>
> On Jun 26, 2013, at 3:18 PM, mar leo wrote:
>
> > Hello!
> > I have a problem with the read_ic.
> >
> > I need to initialize a new vectorial variable for every particle (some kind of normal vector, coming from geometrical issues).. I am currently trying to add a new vectorial block to the Gadget2 read_ic.
> > I read the previous posts about IC and the documentation: from what I understood, the input/output form for Gadget it is the same and they are "managed" in the same "code locations".
> >
> > The fortran routine with which I get the IC file seems to be ok, I assigned the new vector components to every particles like I assigned their position. I can compile successfully, also with the printing I get what I want.
> >
> > The problem comes out when I try to add the new block in Gadget2 (read_ic); in the recent past I have add some scalar quantity (which I would like to have as output), and I had no problem at all. I donīt understand why I have this difficulties with a vector.
> >
> > I updated the following things:
> >
> > allsvars.h
> >
> > -#define IO_NBLOCKS 21
> >
> >
> > -enum iofields /*!< this enumeration lists the defined output blocks in snapshot files. Not all of them need to be present. */
> > {
> > IO_POS,
> > IO_VEL,
> > IO_ID,
> > IO_MASS,
> > IO_U,
> > IO_RHO,
> > .....
> > }; ----> I have 21 fields, as blocks number
> >
> > -Pnorm is defined as:
> > double Pnorm[3] in extern struct particle_data
> >
> > io.c
> >
> > I add the same cases that the variable POS in the various function of io.c (because they are supposed to be similar). See following:
> > case IO_Pnorm:
> > strncpy(Tab_IO_Labels[IO_Pnorm], "Pnorm ", 4);
> > break;
> >
> >
> > I add also the new vectorial quantity in int blockpresent(enum iofields blocknr) , because I donīt want to get any output.
> > #ifndef OUTPUTPnorm
> > if(blocknr == IO_Pnorm)
> > return 0;
> > #endif
> >
> >
> >
> > read_ic.c
> >
> > case IO_Pnorm: /* normal vector */
> > for(n = 0; n < pc; n++)
> > for(k = 0; k < 3; k++)
> > P[offset + n].Pnorm[k] = *fp++;
> > break;
> >
> > I compile and it seems to be fine. But then, when I execute, I get this error message:
> >
> >
> > reading file `../ICs/3Dtropfen' on task=0 (contains 4682 particles.)
> > incorrect block-sizes detected!
> > Task=0 blocknr=3 blksize1=56184 blksize2=0
> > task 0: endrun called with an error level of 1889
> >
> >
> > --------------------------------------------------------------------------
> > MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
> > with errorcode 1889.
> >
> > NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
> > You may or may not see output from other processes, depending on
> > exactly when Open MPI kills them.
> >
> >
> > I am really confused, I donīt understand why I get this and I guess I have not clear what blksize1 and blksize2 actually are...How can I go over the problem?
>
> blksize1 and blksize2 give the size in bytes of the corresponding block of variables that is bracketed by these 4-byte values in the file. In particular, in a file with a correct structure, the values of blksize1 and blksize2 must be equal. Fortran's unformatted I/O produces these block brackets automatically, and gadget uses them too to check that the binary file is correctly read in (or has the right structure to be begin with). Your error could either be caused by producing an incorrectly structured input file, or by modifying the reading routine not in a correct way. In any case, something does't quite match up.
>
> Volker
>
>
>
>
>
>
> > If someone has an idea about it and can explain me why I get this problem only with vectors, would be fantastic.
> > I thanks you a lot in advance.
> >
> > Kind regards,
> >
> > Marzia
> >
> > -----------------------------------------------------------
> >
> > If you wish to unsubscribe from this mailing, send mail to
> > minimalist_at_MPA-Garching.MPG.de with a subject of: unsubscribe gadget-list
> > A web-archive of this mailing list is available here:
> > http://www.mpa-garching.mpg.de/gadget/gadget-list
>
>
>
>
> -----------------------------------------------------------
>
> If you wish to unsubscribe from this mailing, send mail to
> minimalist_at_MPA-Garching.MPG.de with a subject of: unsubscribe gadget-list
> A web-archive of this mailing list is available here:
> http://www.mpa-garching.mpg.de/gadget/gadget-list
Hallo,
first of all thank you very much for your reply, it was helpful. I checked again my fortran routine and it came out that there was, of course, an error I have not seen. I am now able to get a IC that seems to be ok; in practice, I added a last new block after internal energy.
Here you find how it looks like:
---------------------------------------------
Particle Pnorm Typ 0
---------------------------------------------
1 0.0000000 -1.0000000 0.0000000
2 0.0000000 -1.0000000 0.0000000
3 0.0000000 -1.0000000 0.0000000
4 0.0000000 -1.0000000 0.0000000
5 0.0000000 -1.0000000 0.0000000
6 0.0000000 -1.0000000 0.0000000
7 0.0000000 -1.0000000 0.0000000
8 0.0000000 -1.0000000 0.0000000
9 0.0000000 -1.0000000 0.0000000
10 0.0000000 -1.0000000 0.0000000
11 0.0000000 -1.0000000 0.0000000
12 0.0000000 -1.0000000 0.0000000
13 0.0000000 -1.0000000 0.0000000
14 0.0000000 -1.0000000 0.0000000
15 0.0000000 -1.0000000 0.0000000
16 0.0000000 -1.0000000 0.0000000
17 0.0000000 -1.0000000 0.0000000
18 0.0000000 -1.0000000 0.0000000
19 0.0000000 -1.0000000 0.0000000
20 0.0000000 -1.0000000 0.0000000
21 0.0000000 -1.0000000 0.0000000
22 0.0000000 -1.0000000 0.0000000
23 0.0000000 -1.0000000 0.0000000
24 0.0000000 -1.0000000 0.0000000
25 0.0000000 -1.0000000 0.0000000
26 0.0000000 -1.0000000 0.0000000
....... ........
4670 0.0000000 -1.0000000 0.0000000
4671 0.0000000 -1.0000000 0.0000000
4672 0.0000000 -1.0000000 0.0000000
4673 0.0000000 -1.0000000 0.0000000
4674 0.0000000 -1.0000000 0.0000000
4675 0.0000000 -1.0000000 0.0000000
4676 0.0000000 -1.0000000 0.0000000
4677 0.0000000 -1.0000000 0.0000000
4678 0.0000000 -1.0000000 0.0000000
4679 0.0000000 -1.0000000 0.0000000
4680 0.0000000 -1.0000000 0.0000000
4681 0.0000000 -1.0000000 0.0000000
4682 0.0000000 -1.0000000 0.0000000
---------------------------------------------
Particle Pnorm Typ 1
---------------------------------------------
---------------------------------------------
Particle Pnorm Typ 2
---------------------------------------------
---------------------------------------------
Particle Pnorm Typ 3
---------------------------------------------
---------------------------------------------
Particle Pnorm Typ 4
---------------------------------------------
---------------------------------------------
Particle Pnorm Typ 5
--------------------------------------------
But I still have problem. If I try to print this P[i].Pnorm[0], P[i].Pnorm[1], P[i].Pnorm[2] in Gadget
(for example in predict, in the cycle over the whole particles), I do
not get the value I assigned to the particles in the IC files but I get
only 0, 0, 0.
So I thought that probably there is an error also in read_ic but I really do not understand where it could be. My empty_read_buffer in the read_ic is like this:
void empty_read_buffer(enum iofields blocknr, int offset, int pc, int type)
{
int n, k;
float *fp;
#ifdef LONGIDS
long long *ip;
#else
int *ip;
#endif
fp = CommBuffer;
ip = CommBuffer;
switch (blocknr)
{
case IO_POS: /* positions */
for(n = 0; n < pc; n++)
for(k = 0; k < 3; k++)
P[offset + n].Pos[k] = *fp++;
for(n = 0; n < pc; n++)
P[offset + n].Type = type; /* initialize type here as well */
break;
case IO_VEL: /* velocities */
for(n = 0; n < pc; n++)
for(k = 0; k < 3; k++)
P[offset + n].Vel[k] = *fp++;
break;
case IO_ID: /* particle ID */
for(n = 0; n < pc; n++)
P[offset + n].ID = *ip++;
break;
case IO_MASS: /* particle mass */
for(n = 0; n < pc; n++)
P[offset + n].Mass = *fp++;
break;
case IO_U: /* temperature */
for(n = 0; n < pc; n++)
SphP[offset + n].Entropy = *fp++;
break;
case IO_Pnorm: /* normal vector */
for(n = 0; n < pc; n++)
for(k = 0; k < 3; k++)
P[offset + n].Pnorm[k] = *fp++;
break;
case IO_RHO: /* density */
for(n = 0; n < pc; n++)
SphP[offset + n].Density = *fp++;
break;
...
I forgot something, it is clear.. If you have any ideas or suggestion, I would be very happy.
Thanks in advance and kind regards,
Marzia
Received on 2013-07-04 12:28:43