Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add xparity constraints #74

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Add xparity constraints #74

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
803f003
Add xparity
ChristopherHogan Sep 13, 2018
62495d4
Add prototype to maxwell.h
ChristopherHogan Sep 13, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/maxwell/maxwell.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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);
Expand Down
175 changes: 153 additions & 22 deletions src/maxwell/maxwell_constraints.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
}

/**************************************************************************/
Expand Down Expand Up @@ -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;

Expand All @@ -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]);
}
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -302,48 +306,175 @@ 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);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ij2 here should flip i, not j.

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)));
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This isn't right. It's tricky, because we are mirror flipping a vector field, but the vector field is stored in a funny "transverse" coordinate system in order to impose a key mathematical constraint.

The "transverse" coordinate system is specially chosen to have nice properties for mirror flipping in the y and z direction: https://github.com/stevengj/mpb/blob/62495d4f9b3402b31383f0cc255cb5a95d033f7e/src/maxwell/maxwell.c#L398-L399

but it may not be so nice for flipping in x.

}
}
}
}
}

/* 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. */

/* return the number of constant (zero-frequency) bands: */
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!");

Expand All @@ -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))
Expand All @@ -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]);
Expand Down