/* 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");
break;
case 'n':
ed.eps_high = atof(optarg);
CHECK(ed.eps_high > 0.0, "index must be positive");
ed.eps_high = ed.eps_high * ed.eps_high;
break;
case 'f':
ed.eps_high_x = atof(optarg);
CHECK(ed.eps_high_x > 0.0, "fill must be positive");
break;
case 'x':
nx = atoi(optarg);
CHECK(nx > 0, "x size must be positive");
break;
case 'y':
ny = atoi(optarg);
CHECK(ny > 0, "y size must be positive");
break;
case 'z':
nz = atoi(optarg);
CHECK(nz > 0, "z size must be positive");
break;
case 'e':
parity = EVEN_Z_PARITY;
break;
case 'm':
parity = ODD_Z_PARITY;
break;
case 't':
target_freq = fabs(atof(optarg));
do_target = 1;
break;
case 'E':
max_err = fabs(atof(optarg));
CHECK(max_err > 0, "maximum error must be positive");
break;
case 'c':
error_tol = fabs(atof(optarg));
break;
case 'g':
mesh_size = atoi(optarg);
CHECK(mesh_size > 0, "mesh size must be positive");
break;
case '1':
stop1 = 1;
break;
case 'p':
which_preconditioner = 1;
break;
case 'v':
verbose = 1;
break;
default:
usage();
exit(EXIT_FAILURE);
}
if (argc != optind) {
usage();
exit(EXIT_FAILURE);
}
}
#endif
#ifdef ENABLE_PROF
stop1 = 1;
#endif
mesh[0] = mesh[1] = mesh[2] = mesh_size;
printf("Creating Maxwell data...\n");
mdata = create_maxwell_data(nx, ny, nz, &local_N, &N_start, &alloc_N,
num_bands, NUM_FFT_BANDS);
CHECK(mdata, "NULL mdata");
set_maxwell_data_parity(mdata, parity);
printf("Setting k vector to (%g, %g, %g)...\n",
kvector[0], kvector[1], kvector[2]);
update_maxwell_data_k(mdata, kvector, G[0], G[1], G[2]);
printf("Initializing dielectric...\n");
/* set up dielectric function (a simple Bragg mirror) */
set_maxwell_dielectric(mdata, mesh, R, G, epsilon, &ed);
if (verbose && ny == 1 && nz == 1) {
printf("dielectric function:\n");
for (i = 0; i < nx; ++i) {
if (mdata->eps_inv[i].m00 == mdata->eps_inv[i].m11)
printf(" eps(%g) = %g\n", i * 1.0 / nx,
1.0/mdata->eps_inv[i].m00);
else
printf(" eps(%g) = x: %g OR y: %g\n", i * 1.0 / nx,
1.0/mdata->eps_inv[i].m00,
1.0/mdata->eps_inv[i].m11);
}
printf("\n");
}
printf("Allocating fields...\n");
H = create_evectmatrix(nx * ny * nz, 2, num_bands,
local_N, N_start, alloc_N);
Hstart = create_evectmatrix(nx * ny * nz, 2, num_bands,
local_N, N_start, alloc_N);
for (i = 0; i < NWORK; ++i)
W[i] = create_evectmatrix(nx * ny * nz, 2, num_bands,
local_N, N_start, alloc_N);
CHK_MALLOC(eigvals, real, num_bands);
for (iters = 0; iters < PROF_ITERS; ++iters) {
printf("Initializing fields...\n");
for (i = 0; i < H.n * H.p; ++i)
ASSIGN_REAL(Hstart.data[i], rand() * 1.0 / RAND_MAX);
/*****************************************/
if (do_target) {
printf("\nSolving for eigenvectors close to %f...\n", target_freq);
mtdata = create_maxwell_target_data(mdata, target_freq);
op = maxwell_target_operator;
if (which_preconditioner == 1)
pre_op = maxwell_target_preconditioner;
else
pre_op = maxwell_target_preconditioner2;
op_data = (void *) mtdata;
pre_op_data = (void *) mtdata;
}
else {
op = maxwell_operator;
if (which_preconditioner == 1)
pre_op = maxwell_preconditioner;
else
pre_op = maxwell_preconditioner2;
op_data = (void *) mdata;
pre_op_data = (void *) mdata;
}
/*****************************************/
printf("\nSolving for eigenvectors with preconditioning...\n");
evectmatrix_copy(H, Hstart);
eigensolver(H, eigvals,
op, op_data,
pre_op, pre_op_data,
maxwell_parity_constraint, (void *) mdata,
W, NWORK, error_tol, &num_iters, EIGS_DEFAULT_FLAGS);
if (do_target)
eigensolver_get_eigenvals(H, eigvals, maxwell_operator, mdata,
W[0], W[1]);
printf("Solved for eigenvectors after %d iterations.\n", num_iters);
printf("%15s%15s%15s%15s\n","eigenval", "frequency", "exact freq.",
"error");
for (i = 0; i < num_bands; ++i) {
double err;
real freq = sqrt(eigvals[i]);
real exact_freq = bragg_omega(freq, kvector[0], sqrt(ed.eps_high),
ed.eps_high_x, sqrt(ed.eps_low),
1.0 - ed.eps_high_x, 1.0e-7);
printf("%15f%15f%15f%15e\n", eigvals[i], freq, exact_freq,
err = fabs(freq - exact_freq) / exact_freq);
CHECK(err <= max_err, "error exceeds tolerance");
}
printf("\n");
}
if (!stop1) {
/*****************************************/
printf("\nSolving for eigenvectors without preconditioning...\n");
evectmatrix_copy(H, Hstart);
eigensolver(H, eigvals,
op, op_data,
NULL, NULL,
maxwell_parity_constraint, (void *) mdata,
W, NWORK, error_tol, &num_iters, EIGS_DEFAULT_FLAGS);
if (do_target)
eigensolver_get_eigenvals(H, eigvals, maxwell_operator, mdata,
W[0], W[1]);
printf("Solved for eigenvectors after %d iterations.\n", num_iters);
printf("%15s%15s%15s%15s\n","eigenval", "frequency", "exact freq.",
"error");
for (i = 0; i < num_bands; ++i) {
double err;
real freq = sqrt(eigvals[i]);
real exact_freq = bragg_omega(freq, kvector[0], sqrt(ed.eps_high),
ed.eps_high_x, sqrt(ed.eps_low),
1.0 - ed.eps_high_x, 1.0e-7);
printf("%15f%15f%15f%15e\n", eigvals[i], freq, exact_freq,
err = fabs(freq - exact_freq) / exact_freq);
CHECK(err <= max_err, "error exceeds tolerance");
}
printf("\n");
/*****************************************/
printf("\nSolving for eigenvectors without conj. grad...\n");
evectmatrix_copy(H, Hstart);
eigensolver(H, eigvals,
op, op_data,
pre_op, pre_op_data,
maxwell_parity_constraint, (void *) mdata,
W, NWORK - 1, error_tol, &num_iters, EIGS_DEFAULT_FLAGS);
if (do_target)
eigensolver_get_eigenvals(H, eigvals, maxwell_operator, mdata,
W[0], W[1]);
printf("Solved for eigenvectors after %d iterations.\n", num_iters);
printf("%15s%15s%15s%15s\n","eigenval", "frequency", "exact freq.",
"error");
for (i = 0; i < num_bands; ++i) {
double err;
real freq = sqrt(eigvals[i]);
real exact_freq = bragg_omega(freq, kvector[0], sqrt(ed.eps_high),
ed.eps_high_x, sqrt(ed.eps_low),
1.0 - ed.eps_high_x, 1.0e-7);
printf("%15f%15f%15f%15e\n", eigvals[i], freq, exact_freq,
err = fabs(freq - exact_freq) / exact_freq);
CHECK(err <= max_err, "error exceeds tolerance");
}
printf("\n");
/*****************************************/
printf("\nSolving for eigenvectors without precond. or conj. grad...\n");
evectmatrix_copy(H, Hstart);
eigensolver(H, eigvals,
op, op_data,
NULL, NULL,
maxwell_parity_constraint, (void *) mdata,
W, NWORK - 1, error_tol, &num_iters, EIGS_DEFAULT_FLAGS);
if (do_target)
eigensolver_get_eigenvals(H, eigvals, maxwell_operator, mdata,
W[0], W[1]);
printf("Solved for eigenvectors after %d iterations.\n", num_iters);
printf("%15s%15s%15s%15s\n","eigenval", "frequency", "exact freq.",
"error");
for (i = 0; i < num_bands; ++i) {
double err;
real freq = sqrt(eigvals[i]);
real exact_freq = bragg_omega(freq, kvector[0], sqrt(ed.eps_high),
ed.eps_high_x, sqrt(ed.eps_low),
1.0 - ed.eps_high_x, 1.0e-7);
printf("%15f%15f%15f%15e\n", eigvals[i], freq, exact_freq,
err = fabs(freq - exact_freq) / exact_freq);
CHECK(err <= max_err, "error exceeds tolerance");
}
printf("\n");
/*****************************************/
}
destroy_evectmatrix(H);
destroy_evectmatrix(Hstart);
for (i = 0; i < NWORK; ++i)
destroy_evectmatrix(W[i]);
destroy_maxwell_target_data(mtdata);
destroy_maxwell_data(mdata);
free(eigvals);
debug_check_memory_leaks();
return EXIT_SUCCESS;
}
syntax highlighted by Code2HTML, v. 0.9.1