Skip to content

Commit

Permalink
support susceptibilities in user-defined materials (#203)
Browse files Browse the repository at this point in the history 
* trying mkdocs stuff on master branch as readthedocs seems to have trouble with other branches

* synced with stevengj/master

* sync with stevenj master

* updates

* adds support for susceptibilities in user-defined materials; fixes github issue #197

* added --with-susceptibility option to libmeepgeom/user-defined-material test
  • Loading branch information
@stevengj
Homer Reid authored and stevengj committed Feb 16, 2018
1 parent 72705e3 commit 0f6a04b
Show file tree
Hide file tree
Showing 2 changed files with 107 additions and 23 deletions.
35 changes: 26 additions & 9 deletions libmeepgeom/meepgeom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,25 @@ static bool medium_struct_equal(const medium_struct *m1, const medium_struct *m2
susceptibility_list_equal(m1->H_susceptibilities, m2->H_susceptibilities));
}

// garbage collection for susceptibility_list structures.
// Assumes that the 'items' field, if non-empty, was allocated using new[];
// this is automatically the case for python code but is not checked
// for c++ code and will yield runtime errors if a user's user_material_func
// uses e.g. malloc() instead.
static void susceptibility_list_gc(susceptibility_list *sl)
{ if ( !sl || !(sl->num_items) ) return;
delete[] sl->items;
sl->num_items=0;
}

// garbage collection for material structures: called to deallocate memory
// allocated for susceptibilities in user-defined materials.
static void material_gc(material_type m)
{ if ( !m || m->which_subclass!=material_data::MATERIAL_USER ) return;
susceptibility_list_gc( &(m->medium.E_susceptibilities) );
susceptibility_list_gc( &(m->medium.H_susceptibilities) );
}

static bool material_type_equal(const material_type m1, const material_type m2)
{
if (m1 == m2) return true;
Expand Down Expand Up @@ -636,15 +655,7 @@ void geom_epsilon::get_material_pt(material_type &material, const meep::vec &r)
case material_data::MATERIAL_USER:
md->medium = vacuum_medium;
md->user_func(p, md->user_data, &(md->medium));
// TODO: update this to allow user's function to set
// position-dependent susceptibilities. For now
// it's an error if the user's function creates
// any.
if ( (md->medium.E_susceptibilities.num_items>0)
|| (md->medium.H_susceptibilities.num_items>0)
)
meep::abort("susceptibilities in user-defined-materials not yet supported");
return;
return;

// position-independent material or metal: there is nothing to do
case material_data::MEDIUM:
Expand Down Expand Up @@ -675,6 +686,7 @@ double geom_epsilon::chi1p1(meep::field_type ft, const meep::vec &r)
material_type material;
get_material_pt(material, r);
material_epsmu(ft, material, &chi1p1, &chi1p1_inv);
material_gc(material);

return (chi1p1.m00 + chi1p1.m11 + chi1p1.m22)/3;
}
Expand Down Expand Up @@ -829,6 +841,7 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symmetric_matrix *chi1i
get_material_pt(mat, v.center());
trivial:
material_epsmu(meep::type(c), mat, &meps, chi1inv_matrix);
material_gc(mat);
return;
}

Expand Down Expand Up @@ -994,6 +1007,7 @@ void geom_epsilon::fallback_chi1inv_row(meep::component c,
material_type material;
get_material_pt(material, v.center());
material_epsmu(meep::type(c), material, &chi1p1, &chi1p1_inv);
material_gc(material);
if (chi1p1.m01 != 0 || chi1p1.m02 != 0 || chi1p1.m12 != 0
|| chi1p1.m00 != chi1p1.m11 || chi1p1.m11 != chi1p1.m22 ||
chi1p1.m00 != chi1p1.m22) {
Expand Down Expand Up @@ -1138,6 +1152,7 @@ double geom_epsilon::chi(meep::component c, const meep::vec &r, int p) {
chi_val = 0;
};

material_gc(material);
return chi_val;
}

Expand Down Expand Up @@ -1229,6 +1244,7 @@ double geom_epsilon::conductivity(meep::component c, const meep::vec &r) {
default:
cond_val = 0;
}
material_gc(material);

// if the user specified scalar absorbing layers, add their conductivities
// to cond_val (isotropically, for both magnetic and electric conductivity).
Expand Down Expand Up @@ -1315,6 +1331,7 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3],
}
break;
}
material_gc(mat);
}

// add a polarization to the list if it is not already there
Expand Down
95 changes: 81 additions & 14 deletions libmeepgeom/user-defined-material.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,28 @@
#define DATADIR "./"
#endif

/***************************************************************/
/* geometric parameters ****************************************/
/***************************************************************/
#define R1X 0.5 // ellipsoid minor radius
#define R1Y 1.0 // ellipsoid major radius
#define R2 3.0 // cylinder radius

/***************************************************************/
/* parameters in (two-oscillator) Drude-Lorenz model of Ag */
/***************************************************************/
#define EV_UM (1.0/1.23984193)
#define AG_WP (9.01 * EV_UM)
#define AG_F0 0.845

#define AG_FRQ0 1.0e-10
#define AG_GAM0 0.48*EV_UM
#define AG_SIG0 (AG_F0*AG_WP*AG_WP)/(AG_FRQ0*AG_FRQ0)

#define AG_FRQ1 0.065
#define AG_GAM1 0.816*EV_UM
#define AG_SIG1 (AG_F0*AG_WP*AG_WP)/(AG_FRQ1*AG_FRQ1)

using namespace meep;

typedef std::complex<double> cdouble;
Expand Down Expand Up @@ -65,7 +87,6 @@ bool compare_hdf5_datasets(const char *file1, const char *name1,
normDelta = sqrt(normDelta) / size;

double norm12 = fmax( sqrt(norm1), sqrt(norm2) );
printf("Norm{1,2,d} = {%e,%e,%e} (%e)\n",norm1,norm2,normDelta,normDelta/norm12);
if (normDelta > rel_tol*norm12)
return false;
return true;
Expand All @@ -82,25 +103,62 @@ double dummy_eps(const vec &) { return 1.0; }
/* material function that recreates the ellipsoid-in-cylinder */
/* configuration of the cyl-ellipsoid sample code */
/***************************************************************/
typedef struct my_material_func_data
{
double rxInner, ryInner, rOuter;
bool with_susceptibility;
} my_material_func_data;

void my_material_func(vector3 p, void *user_data, meep_geom::medium_struct *m)
{
(void) user_data;
#define R1X 0.5
#define R1Y 1.0
#define R2 3.0
my_material_func_data *data=(my_material_func_data *)user_data;
double rxInner = data->rxInner, rxInner2=rxInner*rxInner;
double ryInner = data->ryInner, ryInner2=ryInner*ryInner;
double rOuter = data->rOuter, rOuter2 =rOuter*rOuter;

double x=p.x, y=p.y;
double x=p.x, x2=x*x, y=p.y, y2=y*y;

// test for point inside inner ellipsoid
double nn;
if ( (x*x/(R1X*R1X) + y*y/(R1Y*R1Y)) < 1.0 )
nn=1.0;
else if ( ( x*x/(R2*R2) + y*y/(R2*R2) ) < 1.0 )
nn=3.5;
else
nn=1.0;
bool innermost = (( x2/rxInner2 + y2/ryInner2 ) < 1.0 );
bool outermost = ( (x*x + y*y) > rOuter2 );
bool in_middle = (!innermost && !outermost);

// set permittivity
double nn = in_middle ? 3.5 : 1.0;
m->epsilon_diag.x = m->epsilon_diag.y = m->epsilon_diag.z = nn*nn;

// add susceptibilities (two-oscillator model for Ag)
if (in_middle&& data->with_susceptibility)
{
m->E_susceptibilities.num_items = 2;
m->E_susceptibilities.items = new meep_geom::susceptibility[2];

m->E_susceptibilities.items[0].sigma_offdiag.x = 0.0;
m->E_susceptibilities.items[0].sigma_offdiag.y = 0.0;
m->E_susceptibilities.items[0].sigma_offdiag.z = 0.0;
m->E_susceptibilities.items[0].sigma_diag.x = AG_SIG0;
m->E_susceptibilities.items[0].sigma_diag.y = AG_SIG0;
m->E_susceptibilities.items[0].sigma_diag.z = AG_SIG0;
m->E_susceptibilities.items[0].frequency = AG_FRQ0;
m->E_susceptibilities.items[0].gamma = AG_GAM0;
m->E_susceptibilities.items[0].noise_amp = 0.0;
m->E_susceptibilities.items[0].drude = true;
m->E_susceptibilities.items[0].is_file = false;

m->E_susceptibilities.items[1].sigma_offdiag.x = 0.0;
m->E_susceptibilities.items[1].sigma_offdiag.y = 0.0;
m->E_susceptibilities.items[1].sigma_offdiag.z = 0.0;
m->E_susceptibilities.items[1].sigma_diag.x = AG_SIG1;
m->E_susceptibilities.items[1].sigma_diag.y = AG_SIG1;
m->E_susceptibilities.items[1].sigma_diag.z = AG_SIG1;
m->E_susceptibilities.items[1].frequency = AG_FRQ1;
m->E_susceptibilities.items[1].gamma = AG_GAM1;
m->E_susceptibilities.items[1].noise_amp = 0.0;
m->E_susceptibilities.items[1].drude = true;
m->E_susceptibilities.items[1].is_file = false;

}

}

/***************************************************************/
Expand All @@ -113,13 +171,16 @@ int main(int argc, char *argv[])
// simple argument parsing
meep::component src_cmpt=Ez;
std::string eps_ref_file = "cyl-ellipsoid-eps-ref.h5";
bool with_susceptibility=true;
for(int narg=1; narg<argc; narg++)
{
if (argv[narg] && !strcmp(argv[narg],"--eps_ref_file"))
{ if (narg+1 == argc)
abort("no option specified for --eps_ref_file");
eps_ref_file=argv[++narg];
}
else if (argv[narg] && !strcmp(argv[narg],"--without_susceptibility"))
with_susceptibility=false;
else
abort("unrecognized command-line option %s",argv[narg]);
};
Expand All @@ -136,8 +197,14 @@ int main(int argc, char *argv[])
: -mirror(X,gv) - mirror(Y,gv);
structure the_structure(gv, dummy_eps, pml(1.0), sym);

my_material_func_data data;
data.with_susceptibility=with_susceptibility;
data.rxInner = R1X;
data.ryInner = R1Y;
data.rOuter = R2;

meep_geom::material_type my_material
= meep_geom::make_user_material(my_material_func, 0);
= meep_geom::make_user_material(my_material_func, (void *)&data);

geometric_object_list g={0, 0};
bool use_anisotropic_averaging=true;
Expand Down

0 comments on commit 0f6a04b

Please sign in to comment.