h5test.cpp

来自mit的fdtd开放性源代码meep

#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <meep.hpp>
#include "meep_internals.hpp"
#include "config.h"
using namespace meep;

const double xsize = 2.0;
const double ysize = 2.0;
const double zsize = 0.6;

const double r = 0.5;
const double eps_k = 2*pi / 1.0;

double funky_eps_2d(const vec &p_) {
  vec p = p_ - vec(xsize / 2, ysize / 2);
  if (fabs(p & p) < r * r)
    return 1.0;
  return 2.0 + cos(p.x() * eps_k) * cos(p.y() * eps_k);
}

double funky_eps_3d(const vec &p_) {
  vec p = p_ - vec(xsize / 2, ysize / 2, zsize / 2);
  if (fabs(p & p) < r * r)
    return 1.0;
  return 2.0 + cos(p.x() * eps_k) * cos(p.y() * eps_k) * cos(p.z() * eps_k);
}

symmetry make_identity(const volume &v)
{
  (void) v; // unused
  return identity();
}

symmetry make_mirrorx(const volume &v)
{
  return mirror(X, v);
}

symmetry make_mirrory(const volume &v)
{
  return mirror(Y, v);
}

symmetry make_mirrorxy(const volume &v)
{
  return mirror(X, v) + mirror(Y, v);
}

symmetry make_rotate4z(const volume &v)
{
  return rotate4(Z, v);
}

typedef symmetry (*symfunc)(const volume &);

const double tol = 1e-8;
double compare(double a, double b, const char *nam, int i0,int i1,int i2) {
  if (fabs(a-b) > 1e-15 + fabs(b) * tol || b != b) {
    master_printf("%g vs. %g differs by\t%g\n", a, b, fabs(a-b));
    master_printf("This gives a fractional error of %g\n", fabs(a-b)/fabs(b));
    abort("Error in %s at (%d,%d,%d)\n", nam, i0,i1,i2);
  }
  return fabs(a-b);
}

double get_reim(complex<double> x, int reim)
{
  return reim ? imag(x) : real(x);
}

bool check_2d(double eps(const vec &), double a, int splitting, symfunc Sf,
	      double kx, double ky,
	      component src_c, int file_c,
	      geometric_volume file_gv,
	      bool real_fields, int expected_rank,
	      const char *name) {
  const volume v = vol2d(xsize, ysize, a);
  structure s(v, eps, no_pml(), Sf(v), splitting);
  fields f(&s);

  f.use_bloch(X, real_fields ? 0.0 : kx);
  f.use_bloch(Y, real_fields ? 0.0 : ky);

  if (real_fields) f.use_real_fields();
  f.add_point_source(src_c, 0.3, 2.0, 0.0, 1.0, v.center(), 1.0, 1);

  if (file_c >= int(Dielectric)) real_fields = true;

  while (f.time() <= 3.0 && !interrupt)
    f.step();

  h5file *file = f.open_h5file(name);
  if (is_derived(file_c))
    f.output_hdf5(derived_component(file_c), file_gv, file);
  else
    f.output_hdf5(component(file_c), file_gv, file);

  file->write("stringtest", "Hello, world!\n");

  delete file;
  all_wait();
  sync();
  file = f.open_h5file(name, h5file::READONLY);

  char *str = file->read("stringtest");
  if (strcmp(str, "Hello, world!\n"))
       abort("Failed to read back string test from %s...", name);

  // compute corner coordinate of file data
  vec loc0(file_gv.get_min_corner());
  ivec iloc0(v.dim);
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    iloc0.set_direction(d, 1+2*int(floor(loc0.in_direction(d)*a-.5)));
    if (file_gv.in_direction(d) == 0.0 &&
	1. - file_gv.in_direction_min(d)*a + 0.5*iloc0.in_direction(d)
	<= 1. + file_gv.in_direction_max(d)*a - 0.5*(iloc0.in_direction(d)+2))
      iloc0.set_direction(d, iloc0.in_direction(d) + 2); // snap to grid
  }
  loc0 = v[iloc0];

  double data_min = infinity, data_max = -infinity;
  double err_max = 0;
  for (int reim = 0; reim < (real_fields ? 1 : 2); ++reim) {
    int rank, dims[2] = {1, 1};

    char dataname[256];
    snprintf(dataname, 256, "%s%s", component_name(file_c),
	     reim ? ".i" : (real_fields ? "" : ".r"));

    double *h5data = file->read(dataname, &rank, dims, 2);
    file->prevent_deadlock(); // hackery
    if (!h5data)
	 abort("failed to read dataset %s:%s\n", name, dataname);
    if (rank != expected_rank)
	 abort("incorrect rank (%d instead of %d) in %s:%s\n",
	       rank, expected_rank, name, dataname);
    if (expected_rank == 1 && 
	file_gv.in_direction_min(X) == file_gv.in_direction_max(X)) {
      dims[1] = dims[0];
      dims[0] = 1;
    }
    vec loc(loc0.dim);
    for (int i0 = 0; i0 < dims[0]; ++i0) {
      for (int i1 = 0; i1 < dims[1]; ++i1) {
	loc.set_direction(X, loc0.in_direction(X) + i0 * v.inva);
	loc.set_direction(Y, loc0.in_direction(Y) + i1 * v.inva);
	int idx = i0 * dims[1] + i1;

	/* Ugh, for rotational symmetries (which mix up components etc.),
	   we can't guarantee that a component is *exactly* the
	   same as its rotated version, and we don't know which one
	   was written to the file. */
	int cs = file_c;
	complex<double> ph = 1.0;
	double diff = fabs(get_reim(f.get_field(file_c, loc), reim) -
			   h5data[idx]);
	for (int sn = 1; sn < f.S.multiplicity(); ++sn) {
	  vec loc2(f.S.transform(loc, sn));
	  int cs2 = f.S.transform(file_c, sn);
	  complex<double> ph2 = f.S.phase_shift(cs2, -sn);
	  double diff2 = fabs(get_reim(f.get_field(cs2, loc2)*ph2, reim) -
			      h5data[idx]);
	  if (diff2 < diff) {
	    loc = loc2;
	    cs = cs2;
	    ph = ph2;
	    diff = diff2;
	  }
	}

	double err = compare(h5data[idx],
			     get_reim(f.get_field(cs, loc) * ph, reim),
			     name, i0,i1,0);
	err_max = max(err, err_max);
	data_min = min(data_min, h5data[idx]);
	data_max = max(data_max, h5data[idx]);
      }
    }
    delete[] h5data;
  }

  file->remove();
  delete file;

  master_printf("Passed %s (%g..%g), err=%g\n", name,
		data_min, data_max,
		err_max / max(fabs(data_min), fabs(data_max)));

  return true;
}

bool check_3d(double eps(const vec &), double a, int splitting, symfunc Sf,
	      component src_c, int file_c,
	      geometric_volume file_gv,
	      bool real_fields, int expected_rank,
	      const char *name) {
  const volume v = vol3d(xsize, ysize, zsize, a);
  structure s(v, eps, no_pml(), Sf(v), splitting);
  fields f(&s);

  if (real_fields) f.use_real_fields();
  f.add_point_source(src_c, 0.3, 2.0, 0.0, 1.0, v.center(), 1.0, 1);

  if (file_c >= Dielectric) real_fields = true;

  while (f.time() <= 3.0 && !interrupt)
    f.step();

  h5file *file = f.open_h5file(name);
  if (is_derived(file_c))
    f.output_hdf5(derived_component(file_c), file_gv, file);
  else
    f.output_hdf5(component(file_c), file_gv, file);

  file->write("stringtest", "Hello, world!\n");

  delete file;
  all_wait();
  sync();
  file = f.open_h5file(name, h5file::READONLY);

  char *str = file->read("stringtest");
  if (strcmp(str, "Hello, world!\n"))
       abort("Failed to read back string test from %s...", name);

  // compute corner coordinate of file data
  vec loc0(file_gv.get_min_corner());
  ivec iloc0(v.dim);
  LOOP_OVER_DIRECTIONS(v.dim, d) {
    iloc0.set_direction(d, 1+2*int(floor(loc0.in_direction(d)*a-.5)));
    if (file_gv.in_direction(d) == 0.0 &&
	1. - file_gv.in_direction_min(d)*a + 0.5*iloc0.in_direction(d)
	<= 1. + file_gv.in_direction_max(d)*a - 0.5*(iloc0.in_direction(d)+2))
      iloc0.set_direction(d, iloc0.in_direction(d) + 2); // snap to grid
  }
  loc0 = v[iloc0];

  double data_min = infinity, data_max = -infinity;
  double err_max = 0;
  for (int reim = 0; reim < (real_fields ? 1 : 2); ++reim) {
    int rank, dims[3] = {1, 1, 1};

    char dataname[256];
    snprintf(dataname, 256, "%s%s", component_name(file_c),
	     reim ? ".i" : (real_fields ? "" : ".r"));