Commit 269b7b0d authored by snuverink_j's avatar snuverink_j
Browse files

Merge branch '469-asciih5block-not-working-with-recent-h5hut-library' into 'master'

Resolve "asciih5block not working with recent H5Hut library?"

Closes #469

See merge request !281
parents 57f4a258 e14d2625
message (STATUS "Configuring BandRF")
set (BANDRF_LIBS
${H5Hut_LIBRARY}
${MPI_CXX_LIBRARIES}
)
add_executable (ascii2h5block ascii2h5block.cpp)
target_link_libraries (ascii2h5block ${BANDRF_LIBS})
install (
TARGETS ascii2h5block
RUNTIME DESTINATION "${CMAKE_INSTALL_PREFIX}/bin"
)
# vi: set et ts=4 sw=4 sts=4:
# Local Variables:
# mode: cmake
# cmake-tab-width: 4
# indent-tabs-mode: nil
# require-final-newline: nil
# End:
\ No newline at end of file
Here a fre comments from Daniel and Hui:
Here a few comments from Daniel and Hui:
Daniel:
Here is a .cpp file I adapted from something that Achim or Jianjun wrote to
generate h5part field files from OPERA 3D field data (*.table text files).
Note that the field data comes without the coordinates, the limits of the
Note that the field data often comes without the coordinates, the limits of the
mesh and the step size in OPERA have to be manually set in this file, and it
is rigged for electric fields or magnetic fields only. So this file is a
total hack, but you should be able to understand what format OPAL needs for
......@@ -15,7 +13,13 @@ the 3D field maps used with BANDRF. If you want electric and magnetic fields
in one h5 file that’s not a problem, just contract the two loops into one
writing both E and H values into the h5 file.
Hui:
My main problem in the past has been that different programs ANSYS, OPERA, COMSOL, ...
all use slightly different default styles for ascii file output,
then you have to recompile ascii2h5block every time, and also for limits of the field,
which are hard-coded in the cpp.
Doing it in python is slower, but not too bad, and removes the need to recompile.
Hui:
Please notice that the 2D mid-plane map is in cylindrical coordinate, while
the 3D map must be in Cartesian coordinate, and be transformed into
......
/*
Purpose: Convert ANSYS E & B-Field data into H5hut (H5block)
format for usage in OPAL.
Usage: ascii2h5block efield.txt hfield.txt ehfout
To visualize use Visit: https://wci.llnl.gov/codes/visit/
Ch. Wang & A. Adelmann, 2011
ToDo: make it more generic
Modification by Chris van Herwaarden and Hui Zhang
The first three rows of a field map that you wish to combine should look like this:
int1 int2 int3
int1 int2 int3
int1 int2 int3
the integers are the amount of steps (or different values) in x y and z
*/
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include "H5hut.h"
int main(int argc,char *argv[]) {
if (argc != 4) {
std::cout << "Wrong number of arguments: ascii2h5block efield.txt (or \"\") hfield.txt (or \"\") ehfout" << std::endl;
//--commlib mpi" << std::endl;
std::exit(1);
}
// // initialize MPI & H5hut
// MPI_Init (&argc, &argv);
// MPI_Comm comm = MPI_COMM_WORLD;
// int comm_size = 1;
// MPI_Comm_size (comm, &comm_size);
// int comm_rank = 0;
// MPI_Comm_rank (comm, &comm_rank);
// H5AbortOnError ();
// H5SetVerbosityLevel (h5_verbosity);
std::string efin(argv[1]);
std::string hfin(argv[2]);
std::string ehfout(argv[3]);
ehfout += std::string(".h5part");
h5_float64_t freq = 72.615 * 1.0e6;
std::cout << "--------------------------------------------------------" << std::endl;
std::cout << "Combine " << efin << " and " << hfin << " to " << ehfout << std::endl;
std::cout << "Frequency " << freq << " [Hz]" << std::endl;
std::ifstream finE, finH;
finH.open(hfin);
finE.open(efin);
bool efield = false, hfield = false;
if (finE.is_open())
efield = true;
else if (efin.empty() == false) {
std::cout << "E-field \"" << efin << "\" could not be opened" << std::endl;
std::exit(1);
}
if (finH.is_open())
hfield = true;
else if (hfin.empty() == false) {
std::cout << "H-field \"" << hfin << "\" could not be opened" << std::endl;
std::exit(1);
}
if (efield == false && hfield == false) {
std::cerr << "Neither E-field \"" << efin
<< "\" nor H-field \"" << hfin
<< "\" could be opened" << std::endl;
std::exit(1);
}
int gridPx = 0, gridPy = 0, gridPz = 0;
char temp[256];
/* Header and grid info */
/* Deletes the first two rows in finE and finH to get rid of header info*/
if (efield) {
finE.getline(temp,256);
finE >> gridPx >> gridPy >> gridPz;
finE.getline(temp,256);
finE.getline(temp,256);
}
int HgridPx = 0, HgridPy = 0, HgridPz = 0;
if (hfield) {
finH.getline(temp,256);
finH >> HgridPx >> HgridPy >> HgridPz;
finH.getline(temp,256);
finH.getline(temp,256);
}
int n = gridPx * gridPy * gridPz; /* number of rows in a column */
int m = HgridPx * HgridPy * HgridPz;
int H5gridPx = std::max(gridPx, HgridPx);
int H5gridPy = std::max(gridPy, HgridPy);
int H5gridPz = std::max(gridPz, HgridPz);
std::cout << "H5block grid" << std::endl;
std::cout << "H5gridPx " << H5gridPx << std::endl;
std::cout << "H5gridPy " << H5gridPy << std::endl;
std::cout << "H5gridPz " << H5gridPz << std::endl;
h5_file_t file = H5OpenFile(ehfout.c_str(), H5_O_WRONLY, H5_PROP_DEFAULT);
if (!file) {
std::cerr << "Could not open output file " << ehfout << std::endl;
std::exit(1);
}
H5SetStep(file, 0);
H5Block3dSetView(file,
0, H5gridPx - 1,
0, H5gridPy - 1,
0, H5gridPz - 1);
std::cout << "number Edata " << n << std::endl;
if (efield) {
h5_float64_t* sEx = new h5_float64_t[n]; /* redefines the sE and sH variables as arrays */
h5_float64_t* sEy = new h5_float64_t[n];
h5_float64_t* sEz = new h5_float64_t[n];
h5_float64_t* FieldstrengthEz = new h5_float64_t[n]; /* redefines the fieldstrength variables as arrays */
h5_float64_t* FieldstrengthEx = new h5_float64_t[n];
h5_float64_t* FieldstrengthEy = new h5_float64_t[n];
double* Ex = new double[n]; /* redefines the E and H variables as arrays */
double* Ey = new double[n];
double* Ez = new double[n];
for (int i = 0; i < n; i++) {
finE >> sEx[i] >> sEy[i] >> sEz[i] >> Ex[i] >> Ey[i] >> Ez[i];
}
finE.close();
h5_float64_t stepEx = (sEx[n-1] - sEx[0]) / (gridPx - 1); /* calculates the stepsizes of the x,y,z of the efield and hfield*/
h5_float64_t stepEy = (sEy[n-1] - sEy[0]) / (gridPy - 1);
h5_float64_t stepEz = (sEz[n-1] - sEz[0]) / (gridPz - 1);
std::cout << "gridPx " << gridPx << " stepEx " << stepEx << std::endl;
std::cout << "gridPy " << gridPy << " stepEy " << stepEy << std::endl;
std::cout << "gridPz " << gridPz << " stepEz " << stepEz << std::endl;
for (int i = 0; i < gridPz; i++) {
for (int j = 0; j < gridPy; j++) {
for (int k = 0; k < gridPx; k++) {
FieldstrengthEx[k + j * gridPx + i * gridPx * gridPy] = static_cast<h5_float64_t>(Ex[i + j * gridPz + k * gridPz * gridPy]);
FieldstrengthEy[k + j * gridPx + i * gridPx * gridPy] = static_cast<h5_float64_t>(Ey[i + j * gridPz + k * gridPz * gridPy]);
FieldstrengthEz[k + j * gridPx + i * gridPx * gridPy] = static_cast<h5_float64_t>(Ez[i + j * gridPz + k * gridPz * gridPy]);
}
}
}
H5Block3dWriteVector3dFieldFloat64 (
file, /*!< IN: file handle */
"Efield", /*!< IN: name of dataset to write */
FieldstrengthEx, /*!< IN: X axis data */
FieldstrengthEy, /*!< IN: Y axis data */
FieldstrengthEz /*!< IN: Z axis data */
);
H5Block3dSetFieldSpacing(file, "Efield", stepEx, stepEy, stepEz);
H5Block3dSetFieldOrigin (file, "Efield", sEx[0], sEy[0], sEz[0]);
}
std::cout << "number Bdata " << m << std::endl;
if (hfield) {
h5_float64_t* sHx = new h5_float64_t[m];
h5_float64_t* sHy = new h5_float64_t[m];
h5_float64_t* sHz = new h5_float64_t[m];
h5_float64_t* FieldstrengthHz = new h5_float64_t[m];
h5_float64_t* FieldstrengthHx = new h5_float64_t[m];
h5_float64_t* FieldstrengthHy = new h5_float64_t[m];
double* Hx = new double[m];
double* Hy = new double[m];
double* Hz = new double[m];
for (int i = 0; i < m; i++) {
finH >> sHx[i] >> sHy[i] >> sHz[i] >> Hx[i] >> Hy[i] >> Hz[i];
}
finH.close();
h5_float64_t stepHx = (sHx[m-1] - sHx[0]) / (HgridPx - 1);
h5_float64_t stepHy = (sHy[m-1] - sHy[0]) / (HgridPy - 1);
h5_float64_t stepHz = (sHz[m-1] - sHz[0]) / (HgridPz - 1);
std::cout << "HgridPx " << HgridPx << " stepHx " << stepHx << std::endl;
std::cout << "HgridPy " << HgridPy << " stepHy " << stepHy << std::endl;
std::cout << "HgridPz " << HgridPz << " stepHz " << stepHz << std::endl;
for (int i = 0; i < HgridPz; i++) {
for (int j = 0; j < HgridPy; j++) {
for (int k = 0; k < HgridPx; k++) {
FieldstrengthHx[k + j * HgridPx + i * HgridPx * HgridPy] = static_cast<h5_float64_t>((Hx[i + j * HgridPz + k * HgridPz * HgridPy]));
FieldstrengthHy[k + j * HgridPx + i * HgridPx * HgridPy] = static_cast<h5_float64_t>((Hy[i + j * HgridPz + k * HgridPz * HgridPy]));
FieldstrengthHz[k + j * HgridPx + i * HgridPx * HgridPy] = static_cast<h5_float64_t>((Hz[i + j * HgridPz + k * HgridPz * HgridPy]));
}
}
}
H5Block3dWriteVector3dFieldFloat64 (
file, /*!< IN: file handle */
"Hfield", /*!< IN: name of dataset to write */
FieldstrengthHx, /*!< IN: X axis data */
FieldstrengthHy, /*!< IN: Y axis data */
FieldstrengthHz /*!< IN: Z axis data */
);
H5Block3dSetFieldSpacing(file, "Hfield", stepHx, stepHy, stepHz);
H5Block3dSetFieldOrigin (file, "Hfield", sHx[0], sHy[0], sHz[0]);
}
H5WriteFileAttribFloat64 (
file, /*!< [in] Handle to open file */
"Resonance Frequency(Hz)", /*!< [in] Name of attribute */
&freq, /*!< [in] Array of attribute values */
1 /*!< [in] Number of array elements */
);
H5CloseFile(file);
std::cout << "Done bye ..." << std::endl;
std::cout << "--------------------------------------------------------" << std::endl;
}
\ No newline at end of file
/*
Purpose: Convert ANSIS E & B-Field data into H5hut (H5block)
format for usage in OPAL.
Usage: ascii2h5block efield.txt hfield.txt ehfield
To visulize use Visit: https://wci.llnl.gov/codes/visit/
Ch. Wang & A. Adelmann, 2011
D. Winklehner, 2013
ToDo: make it more generic / duh -DW
static2: Changed it to reflect new field files from Daniela on Sep. 15 2014
static2a: 0-Hfield for IsoDAR central region RF
*/
#include <fstream>
#include <ios>
#include <iostream>
#include "H5hut.h"
#include <cassert>
#include <string>
#include <algorithm>
#include <cmath>
int main(int argc,char *argv[]) {
if (argc != 3) {
std::cout << "Wrong number of arguments: ascii2h5block efield.txt outfield" << std::endl;
exit(1);
}
std::string efin(argv[1]);
std::string ehfout(argv[2]);
std::string ehfout_c = ehfout + std::string("_CYC.h5part");
std::cout << "--------------------------------------------------------" << std::endl;
std::cout << "Using " << efin << " to create " << ehfout_c << std::endl;
// Open file streams
std::ifstream finE;
finE.open(efin.c_str());
assert(finE.is_open());
// Get number of lines
int nlinesE = std::count(std::istreambuf_iterator<char>(finE),
std::istreambuf_iterator<char>(), '\n');
std::cout << "Lines in finE: " << nlinesE << std::endl;
// Header has 5 lines
int nlines = nlinesE - 5;
// Reset iterator
finE.seekg(0, finE.beg);
std::string templine;
// Skip the 5 header lines
for (int i = 0; i < 5; i++){
std::getline(finE, templine);
}
/* Set frequency (TODO: ask AA if this is the right thing
to do for static fields -DW) */
h5_float64_t freq = 49.2e6; //49.2 MHz ->Hz
h5_float64_t *FieldstrengthEz = new h5_float64_t[nlines];
h5_float64_t *FieldstrengthEx = new h5_float64_t[nlines];
h5_float64_t *FieldstrengthEy = new h5_float64_t[nlines];
h5_float64_t *FieldstrengthHz = new h5_float64_t[nlines];
h5_float64_t *FieldstrengthHx = new h5_float64_t[nlines];
h5_float64_t *FieldstrengthHy = new h5_float64_t[nlines];
// Daniela's latest files don't have the x,y,z coordinates anymore
// Got the following from the readme file accompanying the fields:
h5_float64_t xbegin = -90.0;
h5_float64_t xend = -10.0;
h5_float64_t ybegin = -30.0;
h5_float64_t yend = 30.0;
h5_float64_t zbegin = -1.0;
h5_float64_t zend = 1.0;
// Init the arrays for fields
double *Ex = new double[nlines];
double *Ey = new double[nlines];
double *Ez = new double[nlines];
/*N.B.: Daniela's files now have the structure: Bx,By,Bz; no complex numbers
units are cm, V/cm and Gauss */
for(int i = 0; i < nlines; i++) {
finE >> Ex[i] >> Ey[i] >> Ez[i];
}
finE.close();
double Emax = 0.0;
double E_temp;
int i_temp = 0;
for (int i = 0; i < nlines; i++){
E_temp = std::sqrt(Ex[i] * Ex[i] + Ey[i] * Ey[i] + Ez[i] * Ez[i]);
if (E_temp > Emax) {
Emax = E_temp;
i_temp = i;
}
}
std::cout << "Hardcoded limits: x(" << xbegin << "/" << xend << ") cm" << std::endl;
std::cout << "Hardcoded limits: y(" << ybegin << "/" << yend << ") cm" << std::endl;
std::cout << "Hardcoded limits: z(" << zbegin << "/" << zend << ") cm" << std::endl;
// Set spacing
// TODO: Make program find spacing automatically -DW
double spacing = 0.2;
std::cout << "Hardcoded spacing: " << spacing << " cm" << std::endl;
double gridPx_temp = (xend - xbegin) / spacing + 1.0;
double gridPy_temp = (yend - ybegin) / spacing + 1.0;
double gridPz_temp = (zend - zbegin) / spacing + 1.0;
int gridPx = (int) gridPx_temp;
int gridPy = (int) gridPy_temp;
int gridPz = (int) gridPz_temp;
int nlines_hc = gridPx * gridPy * gridPz;
std::cout << "Hardcoded nlines: " << nlines_hc << std::endl;
std::cout << "File nlines: " << nlines << std::endl;
std::cout << "Grid dimensions: Px = " << gridPx << " , Py = " << gridPy << " , Pz = " << gridPz << std::endl;
std::cout << "E_max = (" << Ex[i_temp] << ", "<< Ey[i_temp] << ", " << Ez[i_temp] << ") V/cm at index " << i_temp << "." << std::endl;
std::cout << "Converting from V/cm and cm to kV/mm and mm before saving h5part" << std::endl;
// Here we also convert from from G to kG
for (int i = 0; i < gridPz; i++) {
for (int j = 0; j < gridPy; j++) {
for (int k = 0; k < gridPx; k++) {
FieldstrengthEx[k+j*gridPx+i*gridPx*gridPy]=0.0; //static_cast<h5_float64_t>(Ex[i+j*gridPz+k*gridPz*gridPy])*1e-4;
FieldstrengthEy[k+j*gridPx+i*gridPx*gridPy]=0.0; //static_cast<h5_float64_t>(Ey[i+j*gridPz+k*gridPz*gridPy])*1e-4;
FieldstrengthEz[k+j*gridPx+i*gridPx*gridPy]=0.0; //static_cast<h5_float64_t>(Ez[i+j*gridPz+k*gridPz*gridPy])*1e-4;
FieldstrengthHx[k+j*gridPx+i*gridPx*gridPy]=static_cast<h5_float64_t>(Ex[i+j*gridPz+k*gridPz*gridPy])*1e-3;
FieldstrengthHy[k+j*gridPx+i*gridPx*gridPy]=static_cast<h5_float64_t>(Ey[i+j*gridPz+k*gridPz*gridPy])*1e-3;
FieldstrengthHz[k+j*gridPx+i*gridPx*gridPy]=static_cast<h5_float64_t>(Ez[i+j*gridPz+k*gridPz*gridPy])*1e-3;
}
}
}
/*
// Here we also convert from V/cm to kV/mm (MV/m)
for (int i = 0; i < gridPz; i++) {
for (int j = 0; j < gridPy; j++) {
for (int k = 0; k < gridPx; k++) {
FieldstrengthEx[k+j*gridPx+i*gridPx*gridPy]=static_cast<h5_float64_t>(Ex[i+j*gridPz+k*gridPz*gridPy])*1e-4;
FieldstrengthEy[k+j*gridPx+i*gridPx*gridPy]=static_cast<h5_float64_t>(Ey[i+j*gridPz+k*gridPz*gridPy])*1e-4;
FieldstrengthEz[k+j*gridPx+i*gridPx*gridPy]=static_cast<h5_float64_t>(Ez[i+j*gridPz+k*gridPz*gridPy])*1e-4;
FieldstrengthHx[k+j*gridPx+i*gridPx*gridPy]=0.0; //static_cast<h5_float64_t>(Hx[i+j*gridPz+k*gridPz*gridPy])*1e-3;
FieldstrengthHy[k+j*gridPx+i*gridPx*gridPy]=0.0; //static_cast<h5_float64_t>(Hy[i+j*gridPz+k*gridPz*gridPy])*1e-3;
FieldstrengthHz[k+j*gridPx+i*gridPx*gridPy]=0.0; //static_cast<h5_float64_t>(Hz[i+j*gridPz+k*gridPz*gridPy])*1e-3;
}
}
}
*/
// Change spacing and limits from cm to mm
spacing *= 10.0;
xbegin *= 10; ybegin *= 10; zbegin *= 10;
xend *= 10; yend *= 10; zend *= 10;
// Write h5part file for OPAL-CYC
h5_err_t h5err;
h5_file_t file = H5OpenFile(ehfout_c.c_str(), H5_O_WRONLY, H5_PROP_DEFAULT);
h5err = H5Block3dSetView(file,
0, gridPx - 1,
0, gridPy - 1,
0, gridPz - 1);
if(file) {
H5SetStep(file, 0);
H5Block3dWriteVector3dFieldFloat64 (
file, /*!< IN: file handle */
"Efield", /*!< IN: name of dataset to write */
FieldstrengthEx, /*!< IN: X axis data */
FieldstrengthEy, /*!< IN: Y axis data */
FieldstrengthEz /*!< IN: Z axis data */
);
h5err = H5Block3dSetFieldSpacing(file, "Efield", spacing, spacing, spacing);
h5err = H5Block3dSetFieldOrigin(file, "Efield", xbegin, ybegin, zbegin);
H5Block3dWriteVector3dFieldFloat64 (
file, /*!< IN: file handle */
"Hfield", /*!< IN: name of dataset to write */
FieldstrengthHx, /*!< IN: X axis data */
FieldstrengthHy, /*!< IN: Y axis data */
FieldstrengthHz /*!< IN: Z axis data */
);
h5err = H5Block3dSetFieldSpacing(file, "Hfield", spacing, spacing, spacing);
h5err = H5Block3dSetFieldOrigin(file, "Hfield", xbegin, ybegin, zbegin);
// Frequency here really in Hz ??? -DW
H5WriteFileAttribFloat64 (
file, /*!< [in] Handle to open file */
"Resonance Frequency(Hz)", /*!< [in] Name of attribute */
&freq, /*!< [in] Array of attribute values */
1 /*!< [in] Number of array elements */
);
H5CloseFile(file);
}
/*
// Write text file for OPAL-T
std::ofstream foutEH;
foutEH.open(ehfout_t.c_str());
assert(foutEH.is_open());
// Change spacing and limits to cm
spacing = spacing * 1.0e-1;
xbegin *= 1.0e-1; ybegin *= 1.0e-1; zbegin *= 1.0e-1;
xend *= 1.0e-1; yend *= 1.0e-1; zend *= 1.0e-1;
// Header
foutEH << "3DDynamic\tXYZ\n";
foutEH << freq*1e-6 << "\n"; // frequency in MHz
foutEH << xbegin << "\t" << xend << "\t" << gridPx-1 << "\n"; //cm / cm / #
foutEH << ybegin << "\t" << yend << "\t" << gridPy-1 << "\n"; //cm / cm / #
foutEH << zbegin << "\t" << zend << "\t" << gridPz-1 << "\n"; //cm / cm / #
// Fielddata
// Here we also convert from V/cm to kV/mm (MV/m) and from G to T (NB the difference to the CYC field!)
foutEH.setf(ios::scientific,ios::floatfield);
foutEH.precision(5);
for (int i=0; i<nlines; i++) {
foutEH << Ex[i]*1e-4 << "\t" << Ey[i]*1e-4 << "\t" << Ez[i]*1e-4 << "\t";
foutEH << Hx[i]*1e-4 << "\t" << Hy[i]*1e-4 << "\t" << Hz[i]*1e-4 << "\n";
}
foutEH.close();
*/
std::cout << "Done, bye ..." << std::endl;
std::cout << "--------------------------------------------------------" << std::endl;
}
......@@ -7,17 +7,22 @@ if (ENABLE_SDDSTOOLS)
endif ()
#add_subdirectory (SDDSReader)
option (ENABLE_MSLANG "Compile MSLang stand-alone compiler" OFF)
option (ENABLE_MSLANG "Compile MSLang stand-alone compiler" OFF)
if (ENABLE_MSLANG)
add_subdirectory (mslang)
endif ()
option (ENABLE_BANDRF "Compile BANDRF field conversion scripts" OFF)
if (ENABLE_BANDRF)
add_subdirectory (BandRF)
endif ()
# vi: set et ts=4 sw=4 sts=4:
# Local Variables:
# mode: cmake
# cmake-tab-width: 4
# indent-tabs-mode: nil
# require-final-newline: nil
# End:
# End:
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment