maxwell_test.c

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

/* Copyright (C) 1999, 2000, 2001, 2002, 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 of the License, 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
 */

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

#include "../src/config.h"
#include <check.h>
#include <blasglue.h>
#include <matrices.h>
#include <eigensolver.h>
#include <maxwell.h>

#if defined(HAVE_GETOPT_H)
#  include <getopt.h>
#endif
#if defined(HAVE_UNISTD_H)
#  include <unistd.h>
#endif

#define NX 32
#define NY 1
#define NZ 1
#define NUM_BANDS 8
#define NUM_FFT_BANDS 5

#define NWORK 3

#define KX 0.5
#define EPS_LOW 1.00
#define EPS_HIGH 9.00
#define EPS_HIGH_X 0.25

#define ERROR_TOL 1e-4

#ifdef ENABLE_PROF
#  define PROF_ITERS 10
#else
#  define PROF_ITERS 1
#endif

#define MESH_SIZE 7

/*************************************************************************/

typedef struct {
     real eps_high, eps_low, eps_high_x;
} epsilon_data;

#define INVERSION_SYM 1

static void epsilon(symmetric_matrix *eps, symmetric_matrix *eps_inv,
		    real r[3], void *edata_v)
{
     epsilon_data *edata = (epsilon_data *) edata_v;
     real eps_val;

#if INVERSION_SYM
     if (fabs(r[0]) < 0.5*edata->eps_high_x 
	 || fabs(r[0]-1.0) < 0.5*edata->eps_high_x)
	  eps_val = edata->eps_high;
#else
     if ((r[0] < edata->eps_high_x && r[0] >= 0.0) ||
	 (r[0] >= 1.0 && r[0] - 1.0 < edata->eps_high_x))
	  eps_val = edata->eps_high;
#endif
     else
	  eps_val = edata->eps_low;
     eps->m00 = eps->m11 = eps->m22 = eps_val;
     eps_inv->m00 = eps_inv->m11 = eps_inv->m22 = 1.0 / eps_val;
#ifdef WITH_HERMITIAN_EPSILON
     CASSIGN_ZERO(eps->m01);
     CASSIGN_ZERO(eps->m02);
     CASSIGN_ZERO(eps->m12);
     CASSIGN_ZERO(eps_inv->m01);
     CASSIGN_ZERO(eps_inv->m02);
     CASSIGN_ZERO(eps_inv->m12);
#else
     eps->m01 = eps->m02 = eps->m12 = 0.0;
     eps_inv->m01 = eps_inv->m02 = eps_inv->m12 = 0.0;
#endif
}

/*************************************************************************/

/* routines for analytic calculation of Bragg bands: */

static const double TWOPI = 6.2831853071795864769252867665590057683943388;

/* We have an analytic expression for k as a function of omega
   for Bragg mirrors.  This will have to be numerically inverted
   to find omega as a function of k.

   n1 and n2 are the indices of the two dielectrics, and f1 and f2
   are their thicknesses as a fraction of the lattice constant
   (we should have f1 + f2 == 1). */
real analytic_bragg_k(real omega, real n1, real f1, real n2, real f2)
{
     real phase1, phase2, c1, s1, c2, s2, b;

     CHECK(fabs(f1 + f2 - 1) < 1e-6, "invalid params to analytic_bragg_k");

     phase1 = TWOPI * n1 * f1 * omega;
     phase2 = TWOPI * n2 * f2 * omega;
     c1 = cos(phase1); s1 = sin(phase1);
     c2 = cos(phase2); s2 = sin(phase2);

     b = c1*c2 - 0.5 * (n1/n2 + n2/n1) * s1*s2;

     if (fabs(b) > 1)
	  return (-1.0);

     return fabs(atan2(sqrt(1-b*b), b) / TWOPI);
}

/* Solve for Bragg omega for the given k and other parameters,
   using omega_guess as a starting guess. 

   We can't do anything clever like Newton's method or even an
   ordinary bisection search because there are regions of omega
   in which analytic_bragg_k is not defined (i.e. in the band gap). */
real bragg_omega(real omega_guess,
		 real k,
		 real n1, real f1, real n2, real f2,
		 real tolerance)
{
     real omega_guess_low = omega_guess - 0.2, 
	  omega_guess_high = omega_guess + 0.2;
     real k_cur;
     real k_best = -1.0, omega_best = 0.0;
     real tol;

     if (omega_guess_low < 0.0)
	  omega_guess_low = 0.0;

     for (tol = (omega_guess_high - omega_guess_low) / 10.0;
	  tol > tolerance;
	  tol *= 0.25) {
	  for (omega_guess = omega_guess_low + tol;
	       omega_guess < omega_guess_high;
	       omega_guess += tol) {
	       k_cur = analytic_bragg_k(omega_guess, n1, f1, n2, f2);
	       if (fabs(k_cur - k) < fabs(k_best - k)) {
		    k_best = k_cur;
		    omega_best = omega_guess;
	       }
	  }

	  CHECK(k_best > 0.0, "No valid omega values in guess range!");

	  omega_guess_low = omega_best - tol;
	  omega_guess_high = omega_best + tol;
     }

     return omega_best;
}


/*************************************************************************/

real norm_diff(scalar *a, scalar *b, int n)
{
     real bmag = 0.0, diffmag = 0.0;
     int i;
     for (i = 0; i < n; ++i) {
	  scalar d;
	  ASSIGN_SCALAR(d,
			SCALAR_RE(b[i]) - SCALAR_RE(a[i]), 
			SCALAR_IM(b[i]) - SCALAR_IM(a[i]));
	  bmag += SCALAR_NORMSQR(b[i]);
	  diffmag += SCALAR_NORMSQR(d);
     }
     return sqrt(diffmag / bmag);
}

/*************************************************************************/

void usage(void)
{
     printf("Syntax: maxwell_test [options]\n"
	    "Options:\n"
            "   -h           Print this help\n"
	    "   -s <seed>    Set random seed\n"
	    "   -k <kx>      Set kx wavevector component [dflt. = %f]\n"
	    "   -b <n>       Compute n bands [default = %d]\n"
	    "   -n <index>   Specify high-dielectric index [dflt. = %f]\n"
	    "   -f <f>       Specify high-index fill fraction [dflt. = %f]\n"
	    "   -x <nx>      Use nx points in x direction [dflt. = %d]\n"
	    "   -y <ny>      Use ny points in y direction [dflt. = %d]\n"
	    "   -z <nz>      Use nz points in z direction [dflt. = %d]\n"
	    "   -e           Solve for TE polarization only.\n"
	    "   -m           Solve for TM polarization only.\n"
	    "   -t <freq>    Set target frequency [dflt. none].\n"
	    "   -c <tol>     Set convergence tolerance [dflt. %e].\n"
	    "   -g <NMESH>   Set mesh size [dflt. %d].\n"
	    "   -1           Stop after first computation.\n"
	    "   -p           Use simple preconditioner.\n"
	    "   -E <err>     Exit with error if the error exceeds <err>\n"
	    "   -v           Verbose output.\n",
	    KX, NUM_BANDS, sqrt(EPS_HIGH), EPS_HIGH_X, NX, NY, NZ,
	    ERROR_TOL, MESH_SIZE);
}

/*************************************************************************/

int main(int argc, char **argv)
{
     maxwell_data *mdata;
     maxwell_target_data *mtdata = NULL;
     int local_N, N_start, alloc_N;
     real R[3][3] = { {1,0,0}, {0,0.01,0}, {0,0,0.01} };
     real G[3][3] = { {1,0,0}, {0,100,0}, {0,0,100} };
     real kvector[3] = {KX,0,0};
     evectmatrix H, Hstart, W[NWORK];
     real *eigvals;
     int i, iters;
     int num_iters;
     int parity = NO_PARITY;
     int nx = NX, ny = NY, nz = NZ;
     int num_bands = NUM_BANDS;
     real target_freq = 0.0;
     int do_target = 0;
     evectoperator op;
     evectpreconditioner pre_op;
     void *op_data, *pre_op_data;
     real error_tol = ERROR_TOL;
     int mesh_size = MESH_SIZE, mesh[3];
     epsilon_data ed;
     int stop1 = 0;
     int verbose = 0;
     int which_preconditioner = 2;
     double max_err = 1e20;

     srand(time(NULL));

     ed.eps_high = EPS_HIGH;
     ed.eps_low = EPS_LOW;
     ed.eps_high_x = EPS_HIGH_X;

#ifdef HAVE_GETOPT
     {
          extern char *optarg;
          extern int optind;
          int c;

          while ((c = getopt(argc, argv, "hs:k:b:n:f:x:y:z:emt:c:g:1pvE:"))
		 != -1)
	       switch (c) {
		   case 'h':
			usage();
			exit(EXIT_SUCCESS);
			break;
		   case 's':
			srand(atoi(optarg));
			break;	
		   case 'k':
			kvector[0] = atof(optarg);
			break;
		   case 'b':
			num_bands = atoi(optarg);
			CHECK(num_bands > 0, "num_bands must be positive");