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

Enable Eigenmode Features with Dispersive Materials #919

Merged
merged 28 commits into from
Jun 28, 2019
Merged

Enable Eigenmode Features with Dispersive Materials #919

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
4582204
first stab
smartalecH May 13, 2019
c08123e
Merge branch 'dispersion' of https://github.com/smartalecH/meep into …
smartalecH May 16, 2019
8ede31f
begin with test files
smartalecH May 22, 2019
b546a9c
bug fixes
smartalecH May 23, 2019
43b6adb
update test
smartalecH Jun 13, 2019
e0ef6ed
merge
smartalecH Jun 14, 2019
23cd5f6
found one bug
smartalecH Jun 15, 2019
a4ac0f2
bug fixes, h5 output, and array-slice output
smartalecH Jun 17, 2019
b0f65c3
fix serial bugs
smartalecH Jun 18, 2019
78d9053
more serial fixes
smartalecH Jun 18, 2019
939a99e
fix mu issue
smartalecH Jun 18, 2019
6155d59
disable h5 support
smartalecH Jun 19, 2019
6eac066
big fix
smartalecH Jun 21, 2019
c4bc10b
generalize matrix inverse and fix parallel bug
smartalecH Jun 21, 2019
97a7b0a
fix memory issue
smartalecH Jun 24, 2019
8a1e635
add tutorial
smartalecH Jun 24, 2019
1f8d91b
Merge branch 'master' of https://github.com/NanoComp/meep into disper…
smartalecH Jun 24, 2019
81f3cc8
minor fixes
smartalecH Jun 24, 2019
041ea9c
refactor
smartalecH Jun 24, 2019
a3e4987
check1
smartalecH Jun 24, 2019
8eebb18
case2
smartalecH Jun 25, 2019
9a33f4c
add docs
smartalecH Jun 25, 2019
2d6aa85
cleanup
smartalecH Jun 25, 2019
6046280
check3
smartalecH Jun 25, 2019
176834a
fix test
smartalecH Jun 25, 2019
6ad454a
1 more fix!
smartalecH Jun 25, 2019
e28c961
updates
smartalecH Jun 25, 2019
8590db1
det fix
smartalecH Jun 26, 2019
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
25 changes: 3 additions & 22 deletions python/geom.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -239,28 +239,9 @@ def transform(self, m):
for s in self.H_susceptibilities:
s.transform(m)

def rotate(self, rotations):
for rot in rotations:
if np.count_nonzero(rot) != 1:
raise ValueError("Each rotation vector should only have 1 coordinate.")
if rot.x != 0: # rotate about x axis
self.transform(Matrix(
Vector3(1,0,0),
Vector3(0,np.cos(rot.x),np.sin(rot.x)),
Vector3(0,-np.sin(rot.x),np.cos(rot.x))
))
elif rot.y != 0: # rotate about z axis
self.transform(Matrix(
Vector3(np.cos(rot.y),0,-np.sin(rot.y)),
Vector3(0,1,0),
Vector3(np.sin(rot.y),0,np.cos(rot.y))
))
else:
self.transform(Matrix(
Vector3(np.cos(rot.z),np.sin(rot.z),0),
Vector3(-np.sin(rot.z),np.cos(rot.z),0),
Vector3(0,0,1)
))
def rotate(self, axis, theta):
T = get_rotation_matrix(axis,theta)
self.transform(T)

def epsilon(self,freq):
return self._get_epsmu(self.epsilon_diag, self.epsilon_offdiag, self.E_susceptibilities, self.D_conductivity_diag, self.D_conductivity_offdiag, freq)
Expand Down
10 changes: 6 additions & 4 deletions python/tests/dispersive_eigenmode.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,10 @@ def call_chi1(self,material,omega):

def verify_output_and_slice(self,material,omega):
# Since the slice routines average the diagonals, we need to do that too:
chi1 = np.square(np.real(np.sqrt(material.epsilon(omega))))
chi1inv = (np.linalg.inv(chi1))
chi1 = material.epsilon(omega).astype(np.complex128)
if np.any(np.imag(chi1) != 0):
chi1 = np.square(np.real(np.sqrt(chi1)))
chi1inv = np.linalg.inv(chi1)
chi1inv = np.diag(chi1inv)
N = chi1inv.size
n = np.sqrt(N/np.sum(chi1inv))
Expand Down Expand Up @@ -101,7 +103,7 @@ def test_chi1_routine(self):
# Now let's rotate LN
import copy
rotLiNbO3 = copy.deepcopy(LiNbO3)
rotLiNbO3.rotate([mp.Vector3(x=np.radians(28)),mp.Vector3(np.radians(45))])
rotLiNbO3.rotate(mp.Vector3(1,1,1),np.radians(34))
self.call_chi1(rotLiNbO3,w0)
self.call_chi1(rotLiNbO3,w1)

Expand Down Expand Up @@ -138,7 +140,7 @@ def test_get_with_dispersion(self):
# Now let's rotate LN
import copy
rotLiNbO3 = copy.deepcopy(LiNbO3)
rotLiNbO3.rotate([mp.Vector3(x=np.radians(28)),mp.Vector3(np.radians(45))])
rotLiNbO3.rotate(mp.Vector3(1,1,1),np.radians(34))
self.verify_output_and_slice(rotLiNbO3,w0)
self.verify_output_and_slice(rotLiNbO3,w1)

Expand Down
2 changes: 1 addition & 1 deletion src/meep.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ const double nan = -7.0415659787563146e103; // ideally, a value never encountere
class h5file;

// Defined in monitor.cpp
void matrix_invert(std::vector<double> &Vinv, std::vector<double> &V);
void matrix_invert(std::complex<double> (&Vinv)[9], std::complex<double> (&V)[9]);

double pml_quadratic_profile(double, void *);

Expand Down
51 changes: 27 additions & 24 deletions src/monitor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -225,29 +225,31 @@ double structure::get_chi1inv(component c, direction d, const ivec &origloc, dou
return 0.0;
}

/* Set Vinv = inverse of V, where both V and Vinv are real-symmetric matrices.*/
void matrix_invert(std::vector<double> &Vinv, std::vector<double> &V) {
/* Set Vinv = inverse of V, where both V and Vinv are complex matrices.*/
void matrix_invert(std::complex<double> (&Vinv)[9], std::complex<double> (&V)[9]) {

double det = (V[0 +3*0] * (V[1 + 3*1]*V[2 +3*2] - V[1 + 3*2]*V[2 +3*1]) -
std::complex<double> det = (V[0 +3*0] * (V[1 + 3*1]*V[2 +3*2] - V[1 + 3*2]*V[2 +3*1]) -
V[0 + 3*1] * (V[0 + 3*1]*V[2 + 3*2] - V[1 + 3*2]*V[0 + 3*2]) +
V[0 + 3*2] * (V[0 + 3*1]*V[1 + 3*2] - V[1 + 3*1]*V[0 + 3*2]));

if (det == 0) abort("meep: Matrix is singular, aborting.\n");

Vinv[0 + 3*0] = 1/det * (V[1 + 3*1]*V[2 + 3*2] - V[1 + 3*2]*V[2 + 3*1]);
Vinv[0 + 3*1] = 1/det * (V[0 + 3*2]*V[2 + 3*1] - V[0 + 3*1]*V[2 + 3*2]);
Vinv[0 + 3*2] = 1/det * (V[0 + 3*1]*V[1 + 3*2] - V[0 + 3*2]*V[1 + 3*1]);
Vinv[1 + 3*0] = 1/det * (V[1 + 3*2]*V[2 + 3*0] - V[1 + 3*0]*V[2 + 3*2]);
Vinv[1 + 3*1] = 1/det * (V[0 + 3*0]*V[2 + 3*2] - V[0 + 3*2]*V[2 + 3*0]);
Vinv[1 + 3*2] = 1/det * (V[0 + 3*2]*V[1 + 3*0] - V[0 + 3*0]*V[1 + 3*2]);
Vinv[2 + 3*0] = 1/det * (V[1 + 3*0]*V[2 + 3*1] - V[1 + 3*1]*V[2 + 3*0]);
Vinv[2 + 3*1] = 1/det * (V[0 + 3*1]*V[2 + 3*0] - V[0 + 3*0]*V[2 + 3*1]);
Vinv[2 + 3*2] = 1/det * (V[0 + 3*0]*V[1 + 3*1] - V[0 + 3*1]*V[1 + 3*0]);
if (det.real() == 0 && det.imag() == 0) abort("meep: Matrix is singular, aborting.\n");
Copy link
Collaborator

Choose a reason for hiding this comment

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

doesn't det == 0.0 work?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I was just being safe -- do you want it changed?

Copy link
Collaborator

Choose a reason for hiding this comment

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

det == 0.0 seems clearer and is safe.

(You can't do det == 0, on the other hand, because C++ doesn't let you mix complex<double> with scalars of other types.)

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Just fixed with 8590db1


Vinv[0 + 3*0] = 1.0/det * (V[1 + 3*1]*V[2 + 3*2] - V[1 + 3*2]*V[2 + 3*1]);
Vinv[0 + 3*1] = 1.0/det * (V[0 + 3*2]*V[2 + 3*1] - V[0 + 3*1]*V[2 + 3*2]);
Vinv[0 + 3*2] = 1.0/det * (V[0 + 3*1]*V[1 + 3*2] - V[0 + 3*2]*V[1 + 3*1]);
Vinv[1 + 3*0] = 1.0/det * (V[1 + 3*2]*V[2 + 3*0] - V[1 + 3*0]*V[2 + 3*2]);
Vinv[1 + 3*1] = 1.0/det * (V[0 + 3*0]*V[2 + 3*2] - V[0 + 3*2]*V[2 + 3*0]);
Vinv[1 + 3*2] = 1.0/det * (V[0 + 3*2]*V[1 + 3*0] - V[0 + 3*0]*V[1 + 3*2]);
Vinv[2 + 3*0] = 1.0/det * (V[1 + 3*0]*V[2 + 3*1] - V[1 + 3*1]*V[2 + 3*0]);
Vinv[2 + 3*1] = 1.0/det * (V[0 + 3*1]*V[2 + 3*0] - V[0 + 3*0]*V[2 + 3*1]);
Vinv[2 + 3*2] = 1.0/det * (V[0 + 3*0]*V[1 + 3*1] - V[0 + 3*1]*V[1 + 3*0]);
}

double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, double omega) const {
smartalecH marked this conversation as resolved.
Show resolved Hide resolved
double res = 0.0;
if (is_mine()){
if (omega == 0)
return chi1inv[c][d] ? chi1inv[c][d][idx] : (d == component_direction(c) ? 1.0 : 0);
// ----------------------------------------------------------------- //
// ---- Step 1: Get instantaneous chi1 tensor ----------------------
// ----------------------------------------------------------------- //
Expand All @@ -268,22 +270,23 @@ double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, dou
my_stuff = B_stuff;
}

std::vector<double> chi1_inv_tensor(9,0);
std::vector<double> chi1_tensor(9,0);
std::complex<double> chi1_inv_tensor[9] = {std::complex<double>(1, 0),std::complex<double>(0, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(1, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(0, 0),std::complex<double>(1, 0)
};
std::complex<double> chi1_tensor[9] = {std::complex<double>(1, 0),std::complex<double>(0, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(1, 0),std::complex<double>(0, 0),
std::complex<double>(0, 0),std::complex<double>(0, 0),std::complex<double>(1, 0)
};

// Set up the chi1inv tensor
// Set up the chi1inv tensor with the DC components
for (int com_it=0; com_it<3;com_it++){
for (int dir_int=0;dir_int<3;dir_int++){
if (chi1inv[comp_list[com_it]][dir_int] )
chi1_inv_tensor[com_it + 3*dir_int] = chi1inv[comp_list[com_it]][dir_int][idx];
else if(dir_int == component_direction(comp_list[com_it]))
chi1_inv_tensor[com_it + 3*dir_int] = 1;
else
chi1_inv_tensor[com_it + 3*dir_int] = 0;
}
}
smartalecH marked this conversation as resolved.
Show resolved Hide resolved


matrix_invert(chi1_tensor, chi1_inv_tensor); // We have the inverse, so let's invert it.

// ----------------------------------------------------------------- //
Expand All @@ -293,7 +296,7 @@ double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, dou
// loop over tensor elements
for (int com_it=0; com_it<3;com_it++){
for (int dir_int=0;dir_int<3;dir_int++){
std::complex<double> eps(chi1_tensor[com_it + 3*dir_int],0.0);
std::complex<double> eps = chi1_tensor[com_it + 3*dir_int];
component cc = comp_list[com_it];
direction dd = (direction)dir_int;
// Loop through and add up susceptibility contributions
Expand Down Expand Up @@ -326,7 +329,7 @@ double structure_chunk::get_chi1inv_at_pt(component c, direction d, int idx, dou
// ----------------------------------------------------------------- //

matrix_invert(chi1_inv_tensor, chi1_tensor); // We have the inverse, so let's invert it.
res = chi1_inv_tensor[component_index(c) + 3*d];
res = chi1_inv_tensor[component_index(c) + 3*d].real();
}
return res;
}
Expand Down