diff --git a/src/maxwell/maxwell.h b/src/maxwell/maxwell.h index 6238d829..d093efd3 100644 --- a/src/maxwell/maxwell.h +++ b/src/maxwell/maxwell.h @@ -80,6 +80,8 @@ typedef struct { #define ODD_Z_PARITY (1<<1) #define EVEN_Y_PARITY (1<<2) #define ODD_Y_PARITY (1<<3) +#define EVEN_X_PARITY (1<<4) +#define ODD_X_PARITY (1<<5) #define MAX_NPLANS 32 @@ -215,6 +217,7 @@ extern void maxwell_ucross_op(evectmatrix Xin, evectmatrix Xout, extern void maxwell_parity_constraint(evectmatrix X, void *data); extern void maxwell_zparity_constraint(evectmatrix X, void *data); extern void maxwell_yparity_constraint(evectmatrix X, void *data); +extern void maxwell_xparity_constraint(evectmatrix X, void *data); extern int maxwell_zero_k_num_const_bands(evectmatrix X, maxwell_data *d); extern void maxwell_zero_k_set_const_bands(evectmatrix X, maxwell_data *d); diff --git a/src/maxwell/maxwell_constraints.c b/src/maxwell/maxwell_constraints.c index a63e6a78..d8c58b91 100644 --- a/src/maxwell/maxwell_constraints.c +++ b/src/maxwell/maxwell_constraints.c @@ -27,7 +27,7 @@ /**************************************************************************/ -/* function to call z and y parity constraints, if necessary */ +/* function to call x, y, and z parity constraints, if necessary */ void maxwell_parity_constraint(evectmatrix X, void *data) { maxwell_data *d = (maxwell_data *) data; @@ -39,6 +39,10 @@ void maxwell_parity_constraint(evectmatrix X, void *data) maxwell_zparity_constraint(X, data); if (d->parity & (EVEN_Y_PARITY | ODD_Y_PARITY)) maxwell_yparity_constraint(X, data); +#ifndef HAVE_MPI + if (d->parity & (EVEN_X_PARITY | ODD_X_PARITY)) + maxwell_xparity_constraint(X, data); +#endif } /**************************************************************************/ @@ -70,7 +74,7 @@ void maxwell_zparity_constraint(evectmatrix X, void *data) int i, j, b, nxy, nz; int zparity = ((d->parity & EVEN_Z_PARITY) ? +1 : ((d->parity & ODD_Z_PARITY) ? -1 : 0)); - + if (zparity == 0) return; @@ -84,12 +88,12 @@ void maxwell_zparity_constraint(evectmatrix X, void *data) else { /* common case (2d system): even/odd == TE/TM */ nxy = d->other_dims * d->last_dim; if (zparity == +1) - for (i = 0; i < nxy; ++i) + for (i = 0; i < nxy; ++i) for (b = 0; b < X.p; ++b) { ASSIGN_ZERO(X.data[(i * X.c + 1) * X.p + b]); } else if (zparity == -1) - for (i = 0; i < nxy; ++i) + for (i = 0; i < nxy; ++i) for (b = 0; b < X.p; ++b) { ASSIGN_ZERO(X.data[(i * X.c) * X.p + b]); } @@ -98,7 +102,7 @@ void maxwell_zparity_constraint(evectmatrix X, void *data) for (i = 0; i < nxy; ++i) { for (j = 0; 2*j <= nz; ++j) { - int ij = i * nz + j; + int ij = i * nz + j; int ij2 = i * nz + (j > 0 ? nz - j : 0); for (b = 0; b < X.p; ++b) { scalar u,v, u2,v2; @@ -156,7 +160,7 @@ double *maxwell_zparity(evectmatrix X, maxwell_data *d) for (i = 0; i < nxy; ++i) for (j = 0; 2*j <= nz; ++j) { - int ij = i * nz + j; + int ij = i * nz + j; int ij2 = i * nz + (j > 0 ? nz - j : 0); for (b = 0; b < X.p; ++b) { scalar u,v, u2,v2; @@ -180,12 +184,12 @@ double *maxwell_zparity(evectmatrix X, maxwell_data *d) mpi_allreduce(zp_scratch, zparity, X.p, double, MPI_DOUBLE, MPI_SUM, mpb_comm); mpi_allreduce(norm_scratch, zp_scratch, X.p, - double, MPI_DOUBLE, MPI_SUM, mpb_comm); + double, MPI_DOUBLE, MPI_SUM, mpb_comm); for (b = 0; b < X.p; ++b) zparity[b] /= zp_scratch[b]; free(zp_scratch); free(norm_scratch); - + return zparity; } @@ -216,7 +220,7 @@ void maxwell_yparity_constraint(evectmatrix X, void *data) for (i = 0; i < nx; ++i) { for (j = 0; 2*j <= ny; ++j) { - int ij = i * ny + j; + int ij = i * ny + j; int ij2 = i * ny + (j > 0 ? ny - j : 0); for (k = 0; k < nz; ++k) { int ijk = ij * nz + k; @@ -273,7 +277,7 @@ double *maxwell_yparity(evectmatrix X, maxwell_data *d) for (i = 0; i < nx; ++i) { for (j = 0; 2*j <= ny; ++j) { - int ij = i * ny + j; + int ij = i * ny + j; int ij2 = i * ny + (j > 0 ? ny - j : 0); for (k = 0; k < nz; ++k) { int ijk = ij * nz + k; @@ -302,20 +306,147 @@ double *maxwell_yparity(evectmatrix X, maxwell_data *d) mpi_allreduce(yp_scratch, yparity, X.p, double, MPI_DOUBLE, MPI_SUM, mpb_comm); mpi_allreduce(norm_scratch, yp_scratch, X.p, - double, MPI_DOUBLE, MPI_SUM, mpb_comm); + double, MPI_DOUBLE, MPI_SUM, mpb_comm); for (b = 0; b < X.p; ++b) yparity[b] /= yp_scratch[b]; free(yp_scratch); free(norm_scratch); - + return yparity; } +/* Similar to the zparity functions above, but for the x -> -x mirror flip. */ + +/* Project X to its even or odd component, so that we can solve + for only one parity of states (the projection operator, like the + mirror flip operator, commutes with the Maxwell operator, so this + projection should not slow convergence). */ +void maxwell_xparity_constraint(evectmatrix X, void *data) +{ +#ifdef HAVE_MPI + CHECK(NULL, "maxwell_xparity_constraint not available when using MPI."); +#endif + maxwell_data *d = (maxwell_data *) data; + int i, j, k, b, nx, ny, nz; + int xparity = ((d->parity & EVEN_X_PARITY) ? +1 : + ((d->parity & ODD_X_PARITY) ? -1 : 0)); + + if (xparity == 0) + return; + + CHECK(d, "null maxwell data pointer!"); + CHECK(X.c == 2, "fields don't have 2 components!"); + + nx = d->nx; + ny = d->ny; + nz = d->nz; + + for (i = 0; 2*i < nx; ++i) { + for (j = 0; j <= ny; ++j) { + int ij = i * ny + j; + int ij2 = i * ny + (j > 0 ? ny - j : 0); + for (k = 0; k < nz; ++k) { + int ijk = ij * nz + k; + int ijk2 = ij2 * nz + k; + for (b = 0; b < X.p; ++b) { + scalar u,v, u2,v2; + u = X.data[(ijk * 2) * X.p + b]; + v = X.data[(ijk * 2 + 1) * X.p + b]; + u2 = X.data[(ijk2 * 2) * X.p + b]; + v2 = X.data[(ijk2 * 2 + 1) * X.p + b]; + ASSIGN_SCALAR(X.data[(ijk * 2) * X.p + b], + 0.5*(SCALAR_RE(u) - xparity*SCALAR_RE(u2)), + 0.5*(SCALAR_IM(u) - xparity*SCALAR_IM(u2))); + ASSIGN_SCALAR(X.data[(ijk * 2 + 1) * X.p + b], + 0.5*(SCALAR_RE(v) + xparity*SCALAR_RE(v2)), + 0.5*(SCALAR_IM(v) + xparity*SCALAR_IM(v2))); + ASSIGN_SCALAR(X.data[(ijk2 * 2) * X.p + b], + 0.5*(SCALAR_RE(u2) - xparity*SCALAR_RE(u)), + 0.5*(SCALAR_IM(u2) - xparity*SCALAR_IM(u))); + ASSIGN_SCALAR(X.data[(ijk2 * 2 + 1) * X.p + b], + 0.5*(SCALAR_RE(v2) + xparity*SCALAR_RE(v)), + 0.5*(SCALAR_IM(v2) + xparity*SCALAR_IM(v))); + } + } + } + } +} + +/* Compute the parity of all of the states in X, returning an array + of the parities (which the caller should deallocate with free). + The parity of an arbitrary state is defined as the expectation value + of the mirror flip operator, and will be +1/-1 for even/odd eigenstates + and something in between for everything else. Assumes that the + columns of X are normalized to 1. */ +double *maxwell_xparity(evectmatrix X, maxwell_data *d) +{ +#ifdef HAVE_MPI + CHECK(NULL, "maxwell_xparity_constraint not available when using MPI."); +#endif + int i, j, k, b, nx, ny, nz; + double *xparity, *xp_scratch, *norm_scratch; + + CHECK(d, "null maxwell data pointer!"); + CHECK(X.c == 2, "fields don't have 2 components!"); + + CHK_MALLOC(xparity, double, X.p); + CHK_MALLOC(xp_scratch, double, X.p); + for (b = 0; b < X.p; ++b) + xp_scratch[b] = 0.0; + CHK_MALLOC(norm_scratch, double, X.p); + for (b = 0; b < X.p; ++b) + norm_scratch[b] = 0.0; + + nx = d->nx; + ny = d->ny; + nz = d->nz; + + for (i = 0; 2*i < nx; ++i) { + for (j = 0; j <= ny; ++j) { + int ij = i * ny + j; + int ij2 = i * ny + (j > 0 ? ny - j : 0); + for (k = 0; k < nz; ++k) { + int ijk = ij * nz + k; + int ijk2 = ij2 * nz + k; + for (b = 0; b < X.p; ++b) { + scalar u,v, u2,v2; + u = X.data[(ijk * 2) * X.p + b]; + v = X.data[(ijk * 2 + 1) * X.p + b]; + u2 = X.data[(ijk2 * 2) * X.p + b]; + v2 = X.data[(ijk2 * 2 + 1) * X.p + b]; + xp_scratch[b] += (ijk == ijk2 ? 1.0 : 2.0) * + (SCALAR_RE(v) * SCALAR_RE(v2) + + SCALAR_IM(v) * SCALAR_IM(v2) - + SCALAR_RE(u) * SCALAR_RE(u2) - + SCALAR_IM(u) * SCALAR_IM(u2)); + norm_scratch[b] += (ijk == ijk2 ? 1.0 : 2.0) * + (SCALAR_RE(v) * SCALAR_RE(v) + + SCALAR_IM(v) * SCALAR_IM(v) + + SCALAR_RE(u) * SCALAR_RE(u) + + SCALAR_IM(u) * SCALAR_IM(u)); + } + } + } + } + + for (b = 0; b < X.p; ++b) + xparity[b] += xp_scratch[b]; + for (b = 0; b < X.p; ++b) + xp_scratch[b] += norm_scratch[b]; + for (b = 0; b < X.p; ++b) + xparity[b] /= xp_scratch[b]; + + free(xp_scratch); + free(norm_scratch); + + return xparity; +} + /**************************************************************************/ /* to fix problems with slow convergence for k ~ 0, manually "put in" the k = 0 solution: first two bands are constant and higher bands are - orthogonal. Note that in the TE/TM case, only one band is constant. + orthogonal. Note that in the TE/TM case, only one band is constant. Also note that, in Fourier space, a constant field corresponds to 1 in the DC component and 0 elsewhere. */ @@ -323,27 +454,27 @@ double *maxwell_yparity(evectmatrix X, maxwell_data *d) int maxwell_zero_k_num_const_bands(evectmatrix X, maxwell_data *d) { int num_const_bands, m_band = 1, n_band = 1; - + CHECK(d, "null maxwell data pointer!"); CHECK(X.c == 2, "fields don't have 2 components!"); - + if (d->parity & (ODD_Z_PARITY | EVEN_Y_PARITY)) m_band = 0; if (d->parity & (ODD_Y_PARITY | EVEN_Z_PARITY)) n_band = 0; - + num_const_bands = m_band + n_band; - + if (num_const_bands > X.p) num_const_bands = X.p; - + return num_const_bands; } void maxwell_zero_k_set_const_bands(evectmatrix X, maxwell_data *d) { int i, j, num_const_bands, m_band = 1, n_band = 1; - + CHECK(d, "null maxwell data pointer!"); CHECK(X.c == 2, "fields don't have 2 components!"); @@ -357,10 +488,10 @@ void maxwell_zero_k_set_const_bands(evectmatrix X, maxwell_data *d) for (j = 0; j < num_const_bands; ++j) { ASSIGN_ZERO(X.data[i * X.p + j]); } - + if (X.Nstart > 0) return; /* DC frequency is not on this process */ - + /* Set DC components to 1 (in two parities) for num_const_bands: */ if (d->parity & (ODD_Z_PARITY | EVEN_Y_PARITY)) @@ -386,7 +517,7 @@ void maxwell_zero_k_constraint(evectmatrix X, void *data) if (X.Nstart > 0) return; /* DC frequency is not on this process */ - + for (j = 0; j < X.p; ++j) { ASSIGN_ZERO(X.data[j]); ASSIGN_ZERO(X.data[X.p + j]);