Dear all,

is there a nice way to create an output for 1 specific grid point?

I want to generate a long term time evolution for a specific value SWTNS::Xi(,j,k) for only one specific grid point (e.g. (i,j,k)=(27,17,27)) without collecting data for an entire line as IOASCII::out1D would do.

Does anyone know how this can be achieved?

Thanks a lot and best regards,

Severin Frank

Thank you Roland!

My simulation is running with the Llama multipatch coordinates (Thornburg04). Your suggested method works fine for points inside the cartesian part!

However if I want to choose a point inside the spherical part (outside of the cartesian grid) I can not gather any data.

I assume that the output thorn is only aware of the cartesian coordinates? Do you (or anyone else) know if there is a way to bypass this issue?

/My current grid setup is:// /

//

/Carpet::domain_from_multipatch            = yes// //CartGrid3D::type                                      = "multipatch"// //CartGrid3D::set_coordinate_ranges_on   = "all maps"// //Coordinates::coordinate_system             = "Thornburg04"// / / //Coordinates::h_cartesian                         = 0.1// //Coordinates::h_radial                               = 0.11// / / //Coordinates::sphere_inner_radius           = 2 // //Coordinates::sphere_outer_radius          = 8 // /

//

/Coordinates::n_angular = 28/

Thanks again and best regards,

Severin Frank

On 15 Oct 2019, at 22:05, Haas, Roland <rhaas@illinois.edumailto:rhaas@illinois.edu> wrote:

Hello Severin,

However if I want to choose a point inside the spherical part (outside of the cartesian grid) I can not gather any data.

I assume that the output thorn is only aware of the cartesian coordinates? Do you (or anyone else) know if there is a way to bypass this issue? The output thorn is CarpetIOASCII, though I am not familiar enough with it to know exactly how it finds out which point to output. It may well not be aware of Llama (it would then use the "local" coordinates instead of the "global", Cartesian ones).

You could write some C code to use Cactus's interpolator to get interpolate the data. If the coordinate given is the location of an actual grid point then there is (effectively) no interpolation. For an example on how to use it, please see the Cactus docs:

https://www.einsteintoolkit.org/referencemanual/ReferenceManualch2.html#x4-1...

You could (ab)use CactusNumerical's InterpToArray thorn which lets you do the same without writing code. For an example use, please see its tests though for your use you will want to use the "scalar_vars" etc parameters instead of "array1d_vars".

Hi,

You can use CarpetIOASCII and CarpetIOHDF5 to do what you want, more or less.

Each of the 6 angular patches will be output in separate files. Each one has a mapping x,y,z->rho,sigma,r, where rho and sigma are two of the angular coordinates (different for each patch), and r is *always* the radial coordinate. So, if you want to output a line in the radial direction, you need to use

CarpetIOASCII::out1d_vars = "myvar" CarpetIOASCII::out1d_z = yes

I believe this is enough. You will get the angular coordinates (0,0) for each patch, corresponding to the +x, -x, +y, -y and +z, -z axes, and the "z" coordinate in the 1D output file will be the radial direction. This *might* depend on the value of n_angular; the angular coordinates chosen will come from out1D_zline_x and out1D_zline_y, which default to 0. If the coordinates are such that 0 is between two grid points, then you might not get any data. This will depend on whether n_angular is even or odd, but I don't remember the details from memory. It's also possible that CarpetIOASCII will give you the nearest point in that case; I'm not sure.

You can do the same trick for output at a constant radius, either 1D or 2D, BUT you have to watch out for the fact that the radius will be chosen as, eg. for 2D output, the value of CarpetIOASCII::out2D_xyplane_z, which defaults to 0, where there is no data since the angular patches start at a finite radius r > 0. So if you want a 2D surface at r = 100, you need to set CarpetIOASCII::out2D_xyplane_z = 100.

I believe the same works with HDF5. It should also work for a single point; 0d output; see the parameters in CarpetIOASCII/param.ccl. Just remember that the radial coordinate will be called "z"!

-- Ian Hinder Research Software Engineer University of Manchester, UK

Hello Erik,

I cannot say that I have much experience with the FileReader.

From the error message you are getting, you seem to be using the PUGH driver, is this correct?

This error is triggered in repos/cactuspughio/IOHDF5Util/src/RecoverVar.c when there is no Cactus variable of that name. Looking at the code in there it seems that it takes the name of the dataset in the file (which is IDBRILLMOL::phi since that is the name of the thorn writing it) an then looks for a Cactus variable of that name. It then (with a warning "Ignoring dataset 'IDBRILLMOL::phi'") skips this dataset.

Looking at the code in IOHDF5Util (and PUGHIOHDF5) it seems to me as if PUGH does not support the "alias" setting (CarpetIOHDF5 does).

To make this work you can either (I think) rename the dataset in the hdf5 file to have the expected name (eg using h5py) or extend RecoverVar.c with something like this (untested!):

--8<-- #ifdef IOUTIL_IOGH_HAS_ALIAS for(vindex = 0 ; vindex < CCTK_NumVars() ; vindex++) { if(ioUtilGH->alias && ioUtilGH->alias[vindex] && CCTK_EQUALS(ioUtilGH->alias[vindex], fullname)) { free(fullname); fullname = strdup(ioUtilGH->alias[vindex]); break; // leaves vindex at correct value } } if(vindex == CCTK_NumVars()) // did not find anything #endif // old code continues below /* check if there is a matching variable */ vindex = CCTK_VarIndex (fullname); --8<--

which more or less mimics what Carpet does in repos/carpet/CarpetIOHDF5/src/Input.cc

Yours, Roland

Thank you!, the method that you give me to stop the evolution, works perfectly. However, there is an error, I can't read the output data as initial data for my Brill wave evolution.

My thorn IDBrillMoL evolves the conformal factor (which I named phi because the initial data is obtained through a wave-like equation), and when the criteria is reached, it stops the evolution and saves a checkpoint with that iteration. For the evolution I wrote a thorn named BrillEvolve and I assigned memory to a variable named "brillpsi" which will store the initial conformal factor. When I read the data (checkpoint) I got the error: No matching variable found for 'IDBRILLMOL::phi'

According to the readingfile interface, I can read a variable and set in a different variable with the option alias, so, I try

IO::filereader_ID_vars = "BrillEvolve::brillpsi{alias='IDBRILLMOL::phi'}"

but I got the error above described. Can you help me with this? Am I misunderstanding how the readingfile interface works?

Erik

El lun., 14 oct. 2019 a las 9:44, Haas, Roland (rhaas@illinois.edu) escribió:

...

Hello Erik,

as to how to decide to stop: the simplest way would be to compute a pointwise stopping criterion (maybe magnitude of the right-hand-side or so) then do a reduction on it using the CCTK_Reduce function (assuming you are using Carpet), see eg. the code in Hydro_Analysis to compute the center of mass:

https://bitbucket.org/einsteintoolkit/einsteinanalysis/src/master/Hydro_Anal...

the function itself is documented here:

https://www.einsteintoolkit.org/usersguide/UsersGuidech9.html#verbatim-63

Once you have decided to stop, terminate the simulation by calling CCTK_TerminateNext ( https://www.einsteintoolkit.org/referencemanual/ReferenceManualch2.html#x4-1... ) making sure to have enabled the option to write a checkpoint when terminating (io::checkpoint_on_terminate = yes).

You can read data from the checkpoint as initial data using the filereader interface

https://www.einsteintoolkit.org/thornguide/CactusBase/IOUtil/documentation.h...

and if you need to change grid structure, you can use the ReadInterpolate thorn (https://github.com/rhaas80/ReadInterpolate/tree/master/ReadInterpolate), which is however a bit more complex to use.

If you need to change data in between what was checkpointed and what needs to be evolved you will have to add routines to the INITIAL bin after IOUtil_RecoverIDFromDatafiles (see IOUtil's schedule.ccl).

Hope this helps getting you started.

Yours, Roland

...

Dear all

I'm trying to construct Brill wave initial data by transforming the elliptic pde equation for the conformal factor in to a wave-like equation with a source. After reaching a stationary state, I want to stop the evolution and give the stationary state as an initial data for the BSSN system.

I'm evolving the wave-like equation with MoL and my question is, How can

I

...

stop the MoL evolution if the solution reaches the numerical tolerance

and

...

then pass it as initial data?

Erik

-- My email is as private as my paper mail. I therefore support encrypting and signing email messages. Get my PGP key from http://pgp.mit.edu .