meep.hpp

麻省理工的计算光子晶体的程序

/* Copyright (C) 2006 Massachusetts Institute of Technology
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2, or (at your option)
%  any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software Foundation,
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#ifndef MEEP_H
#define MEEP_H

#include <stdio.h>
#include <math.h>

#include "meep/vec.hpp"
#include "meep/mympi.hpp"

namespace meep {

extern bool quiet; // if true, suppress all non-error messages from Meep

const double pi = 3.141592653589793238462643383276;

#ifdef INFINITY
const double infinity = INFINITY;
#else
const double infinity = 1e20; // this should be big enough for us
#endif

#ifdef NAN
const double nan = NAN;
#else
const double nan = -1.23454321e123; // ideally, a value never encountered in practice
#endif

class polarizability_identifier {
 public:
  double gamma, omeganot;
  double energy_saturation, saturated_sigma;
  bool operator==(const polarizability_identifier &);
};
class polarizability;
class polarization;
class grace;

// h5file.cpp: HDF5 file I/O.  Most users, if they use this
// class at all, will only use the constructor to open the file, and
// will otherwise use the fields::output_hdf5 functions.
class h5file {
public:
  typedef enum {
    READONLY, READWRITE, WRITE
  } access_mode;
  
  h5file(const char *filename_, access_mode m=READWRITE, bool parallel_=true);
  ~h5file(); // closes the files (and any open dataset)
  
  bool ok();
  
  double *read(const char *dataname, int *rank, int *dims, int maxrank);
  void write(const char *dataname, int rank, const int *dims, double *data,
	     bool single_precision = true);
  
  char *read(const char *dataname);
  void write(const char *dataname, const char *data);
  
  void create_data(const char *dataname, int rank, const int *dims,
		   bool append_data = false,
		   bool single_precision = true);
  void extend_data(const char *dataname, int rank, const int *dims);
  void create_or_extend_data(const char *dataname, int rank,
			     const int *dims,
			     bool append_data, bool single_precision);
  void write_chunk(int rank, const int *chunk_start, const int *chunk_dims,
		   double *data);
  void done_writing_chunks();
  
  void read_size(const char *dataname, int *rank, int *dims, int maxrank);
  void read_chunk(int rank, const int *chunk_start, const int *chunk_dims,
		  double *data);
  
  void remove();
  void remove_data(const char *dataname);
  
  const char *file_name() const { return filename; }

  void prevent_deadlock(); // hackery for exclusive mode
private:
  access_mode mode;
  char *filename;
  bool parallel;

  bool is_cur(const char *dataname);
  void unset_cur();
  void set_cur(const char *dataname, void *data_id);
  char *cur_dataname;
  bool cur_append_data;
  
  /* linked list to keep track of which datasets we are extending...
     this is necessary so that create_or_extend_data can know whether
     to create (overwrite) a dataset or extend it. */
  struct extending_s {
    int dindex;
    char *dataname;
    struct extending_s *next;
  } *extending;
  extending_s *get_extending(const char *dataname) const;
  
  /* store hid_t values as hid_t* cast to void*, so that
     files including meep.h don't need hdf5.h */
  void *id; /* file */
  void *cur_id; /* dataset, if any */

  void *get_id(); // get current (file) id, opening/creating file if needed
  void close_id();
};

#define DEFAULT_SUBPIXEL_TOL 1e-4
#define DEFAULT_SUBPIXEL_MAXEVAL 100000

/* This class is used to compute position-dependent material properties
   like the dielectric function, polarizability sigma, 
   nonlinearities, et cetera.  Simple cases of stateless functions are
   handled by canned subclasses below, but more complicated cases
   can be handled by creating a user-defined subclass of material_function.
   It is useful to group different properties into one class because
   it is likely that complicated implementations will share state between
   properties. */
class material_function {
  material_function(const material_function &ef) {} // prevent copying
public:
  material_function() : omega(nan), gamma(nan), deps(nan), energy_sat(nan) {}
  
  virtual ~material_function() {}
  
  /* Specify a restricted volume: all subsequent eps/sigma/etc
     calls will be for points inside gv, until the next set_volume. */
  virtual void set_volume(const geometric_volume &gv) {}
  virtual void unset_volume(void) {} // unrestrict the volume
  
  /* scalar dielectric function */
  virtual double eps(const vec &r) { return 1.0; }
  
  /* Return interface normal and/or average dielectric eps and 1/eps
     in a given volume gv.   These are virtual so that e.g. libctl
     can override them with more efficient geometry-based routines. */
  virtual vec normal_vector(const geometric_volume &gv);
  virtual void meaneps(double &meps, double &minveps, vec &normal,
		       const geometric_volume &gv, 
		       double tol=DEFAULT_SUBPIXEL_TOL, 
		       int maxeval=DEFAULT_SUBPIXEL_MAXEVAL);
  

  /* polarizability sigma function */
  virtual double sigma(const vec &r) { return 0.0; }
  /* specify polarizability used for subsequent calls to sigma(r) */
  virtual void set_polarizability(double omega_, double gamma_, 
				  double deps_, double energy_sat_) {
    omega = omega_; gamma = gamma_; deps = deps_; energy_sat = energy_sat_;
  }
  
  // Nonlinear susceptibilities
  virtual bool has_chi3() { return false; }
  virtual double chi3(const vec &r) { return 0.0; }
  virtual bool has_chi2() { return false; }
  virtual double chi2(const vec &r) { return 0.0; }
  
  // TODO: dielectric tensor, ...

protected:
  // current polarizability for calls to sigma(r):
  double omega, gamma, deps, energy_sat;
};

class simple_material_function : public material_function {
  double (*f)(const vec &);
  
public:
  simple_material_function(double (*func)(const vec &)) { f = func; }
  
  virtual ~simple_material_function() {}
  
  virtual double eps(const vec &r) { return f(r); }
  virtual double sigma(const vec &r) { return f(r); }
  virtual double chi3(const vec &r) { return f(r); }
  virtual double chi2(const vec &r) { return f(r); }
};

class structure;

class structure_chunk {
 public:
  double a, Courant, dt; // res. a, Courant num., and timestep dt=Courant/a
  double *eps, *chi3[NUM_FIELD_COMPONENTS], *chi2[NUM_FIELD_COMPONENTS];
  double *inveps[NUM_FIELD_COMPONENTS][5];
  double *C[5][NUM_FIELD_COMPONENTS];
  double *Cdecay[5][NUM_FIELD_COMPONENTS][5];
  volume v;  // integer volume that could be bigger than non-overlapping gv below
  geometric_volume gv;
  polarizability *pb;

  int refcount; // reference count of objects using this structure_chunk

  ~structure_chunk();
  structure_chunk(const volume &v,
            const geometric_volume &vol_limit, double Courant, int proc_num);
  structure_chunk(const structure_chunk *);
  void set_epsilon(material_function &eps,
                   bool use_anisotropic_averaging,
		   double tol, int maxeval);
  void set_chi3(material_function &eps);
  void set_chi2(material_function &eps);
  void use_pml(direction, double dx, double boundary_loc, double strength);
  void update_pml_arrays();
  void update_Cdecay();

  void add_polarizability(double sigma(const vec &), double omega, double gamma,
                          double delta_epsilon = 1.0, double energy_saturation = 0.0);
  void add_polarizability(material_function &sigma, double omega, double gamma,
                          double delta_epsilon = 1.0, double energy_saturation = 0.0);

  void mix_with(const structure_chunk *, double);

  int n_proc() const { return the_proc; } // Says which proc owns me!
  int is_mine() const { return the_is_mine; }

  void remove_polarizabilities();

  // monitor.cpp
  double get_inveps(component, direction, const ivec &iloc) const;
  double max_eps() const;
 private:
  double pml_fmin;
  int the_proc;
  int the_is_mine;
};

// linked list of descriptors for boundary regions (currently just for PML)
class boundary_region {
public:
  typedef enum { NOTHING_SPECIAL, PML } boundary_region_kind;
  
  boundary_region() :
    kind(NOTHING_SPECIAL), thickness(0.0), strength(1.0), d(NO_DIRECTION), side(Low), next(0) {}
  boundary_region(boundary_region_kind kind, double thickness, double strength, direction d, boundary_side side, boundary_region *next = 0) : kind(kind), thickness(thickness), strength(strength), d(d), side(side), next(next) {}

  boundary_region(const boundary_region &r) : kind(r.kind), thickness(r.thickness), strength(r.strength), d(r.d), side(r.side) { 
    next = r.next ? new boundary_region(*r.next) : 0;
  }

  ~boundary_region() { if (next) delete next; }
  
  void operator=(const boundary_region &r) {
    kind = r.kind; thickness = r.thickness; strength = r.strength;
    d = r.d; side = r.side;
    if (next) delete next;
    next = r.next ? new boundary_region(*r.next) : 0;
  }
  boundary_region operator+(const boundary_region &r0) const {
    boundary_region r(*this), *cur = &r;
    while (cur->next) cur = cur->next;
    cur->next = new boundary_region(r0);
    return r;
  }

  boundary_region operator*(double strength_mult) const {
    boundary_region r(*this), *cur = &r;
    while (cur) {
      cur->strength *= strength_mult;
      cur = cur->next;
    }
    return r;
  }

  void apply(structure *s) const;
  void apply(const structure *s, structure_chunk *sc) const;

private:
  boundary_region_kind kind;
  double thickness, strength;
  direction d;
  boundary_side side;
  boundary_region *next;
};

boundary_region pml(double thickness, direction d, boundary_side side);
boundary_region pml(double thickness, direction d);
boundary_region pml(double thickness);
#define no_pml() boundary_region()

class structure {
 public:
  structure_chunk **chunks;
  int num_chunks;
  volume v, user_volume;
  double a, Courant, dt; // res. a, Courant num., and timestep dt=Courant/a
  geometric_volume gv;
  symmetry S;
  const char *outdir;
  volume *effort_volumes;
  double *effort;
  int num_effort_volumes;

  ~structure();
  structure();
  structure(const volume &v, material_function &eps,
	    const boundary_region &br = boundary_region(),
	    const symmetry &s = meep::identity(),
	    int num_chunks = 0, double Courant = 0.5,
	    bool use_anisotropic_averaging=false,
	    double tol=DEFAULT_SUBPIXEL_TOL,
	    int maxeval=DEFAULT_SUBPIXEL_MAXEVAL);
  structure(const volume &v, double eps(const vec &), 
	    const boundary_region &br = boundary_region(),
	    const symmetry &s = meep::identity(),
	    int num_chunks = 0, double Courant = 0.5,
	    bool use_anisotropic_averaging=false,
	    double tol=DEFAULT_SUBPIXEL_TOL,
	    int maxeval=DEFAULT_SUBPIXEL_MAXEVAL);
  structure(const structure *);
  structure(const structure &);

  void set_materials(material_function &mat,
		     bool use_anisotropic_averaging=true,
		     double tol=DEFAULT_SUBPIXEL_TOL,
		     int maxeval=DEFAULT_SUBPIXEL_MAXEVAL);
  void set_epsilon(material_function &eps,
                   bool use_anisotropic_averaging=true,
		   double tol=DEFAULT_SUBPIXEL_TOL,
		   int maxeval=DEFAULT_SUBPIXEL_MAXEVAL);
  void set_epsilon(double eps(const vec &),
                   bool use_anisotropic_averaging=true,
		   double tol=DEFAULT_SUBPIXEL_TOL,
		   int maxeval=DEFAULT_SUBPIXEL_MAXEVAL);
  void set_chi3(material_function &eps);
  void set_chi3(double eps(const vec &));
  void set_chi2(material_function &eps);
  void set_chi2(double eps(const vec &));
  polarizability_identifier
     add_polarizability(double sigma(const vec &), double omega, double gamma,
                  double delta_epsilon = 1.0, double energy_saturation = 0.0);
  polarizability_identifier
     add_polarizability(material_function &sigma, double omega, double gamma,
                  double delta_epsilon = 1.0, double energy_saturation = 0.0);

  void remove_polarizabilities();

  void set_output_directory(const char *name);
  void mix_with(const structure *, double);

  bool equal_layout(const structure &) const;
  void print_layout(void) const;

  // monitor.cpp
  double get_inveps(component, direction, const ivec &origloc) const;
  double get_inveps(component, direction, const vec &loc) const;
  double get_eps(const vec &loc) const;
  double max_eps() const;

  friend class boundary_region;

 private: