diff --git a/grudge/function_registry.py b/grudge/function_registry.py index df709f93e7ab13fd180d91cca0cce30735ae8f87..31961b6f48a26183b09c6673e04e860eeb7128d1 100644 --- a/grudge/function_registry.py +++ b/grudge/function_registry.py @@ -88,6 +88,7 @@ class CElementwiseUnaryFunction(Function): i = Variable("i") if self.identifier == "fabs": # FIXME + # Loopy has a type-adaptive "abs", but no "fabs". func_name = "abs" cached_name += func_name diff --git a/grudge/models/em.py b/grudge/models/em.py index bf7495e2c9447b54d37d49a78905e549d9df54d8..3422c3be14f8ec10d872d46adc29caaa155131dd 100644 --- a/grudge/models/em.py +++ b/grudge/models/em.py @@ -29,8 +29,6 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from six.moves import range - from pytools import memoize_method from grudge.models import HyperbolicOperator @@ -110,7 +108,6 @@ class MaxwellOperator(HyperbolicOperator): self.current = current self.incident_bc_data = incident_bc - def flux(self, w): """The template for the numerical flux for variable coefficients. @@ -127,15 +124,14 @@ class MaxwellOperator(HyperbolicOperator): epsilon = self.epsilon mu = self.mu - - Z_int = (mu/epsilon)**0.5 - Y_int = 1/Z_int - Z_ext = (mu/epsilon)**0.5 - Y_ext = 1/Z_ext + Z_int = (mu/epsilon)**0.5 # noqa: N806 + Y_int = 1/Z_int # noqa: N806 + Z_ext = (mu/epsilon)**0.5 # noqa: N806 + Y_ext = 1/Z_ext # noqa: N806 if self.flux_type == "lf": - if self.fixed_material: - max_c = (self.epsilon*self.mu)**(-0.5) + # if self.fixed_material: + # max_c = (self.epsilon*self.mu)**(-0.5) return join_fields( # flux e, @@ -178,14 +174,13 @@ class MaxwellOperator(HyperbolicOperator): e, h = self.split_eh(w) nabla = sym.nabla(self.dimensions) + def e_curl(field): return self.space_cross_e(nabla, field) def h_curl(field): return self.space_cross_h(nabla, field) - - # in conservation form: u_t + A u_x = 0 return join_fields( (self.current - h_curl(h)), @@ -199,7 +194,6 @@ class MaxwellOperator(HyperbolicOperator): pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e)) pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h)) - return join_fields(-pec_e, pec_h) def pmc_bc(self, w): @@ -218,15 +212,14 @@ class MaxwellOperator(HyperbolicOperator): absorb_normal = sym.normal(self.absorb_tag, self.dimensions) - e, h = self.split_eh(w) if self.fixed_material: epsilon = self.epsilon mu = self.mu - absorb_Z = (mu/epsilon)**0.5 - absorb_Y = 1/absorb_Z + absorb_Z = (mu/epsilon)**0.5 # noqa: N806 + absorb_Y = 1/absorb_Z # noqa: N806 absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e)) absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h)) @@ -284,7 +277,6 @@ class MaxwellOperator(HyperbolicOperator): (self.incident_tag, self.incident_bc(w)), ] - def flux(pair): return sym.interp(pair.dd, "all_faces")(self.flux(pair)) @@ -294,10 +286,8 @@ class MaxwellOperator(HyperbolicOperator): flux(sym.int_tpair(w)) + sum( flux(sym.bv_tpair(tag, w, bc)) - for tag, bc in tags_and_bcs) - ) )) / material_divisor - - + for tag, bc in tags_and_bcs) + ))) / material_divisor @memoize_method def partial_to_eh_subsets(self): @@ -313,31 +303,12 @@ class MaxwellOperator(HyperbolicOperator): return tuple(partial_to_all_subset_indices( [e_subset, h_subset])) - def split_eps_mu_eh(self, w): - """Splits an array into epsilon, mu, E and H components. - - Only used for fluxes. - """ - e_idx, h_idx = self.partial_to_eh_subsets() - epsilon, mu, e, h = w[[0]], w[[1]], w[e_idx+2], w[h_idx+2] - - from grudge.flux import FluxVectorPlaceholder as FVP - if isinstance(w, FVP): - return ( - FVP(scalars=epsilon), - FVP(scalars=mu), - FVP(scalars=e), - FVP(scalars=h)) - else: - return epsilon, mu, make_obj_array(e), make_obj_array(h) - def split_eh(self, w): "Splits an array into E and H components" e_idx, h_idx = self.partial_to_eh_subsets() e, h = w[e_idx], w[h_idx] - return e,h - #return make_obj_array(e), make_obj_array(h) + return e, h def get_eh_subset(self): """Return a 6-tuple of :class:`bool` objects indicating whether field @@ -384,8 +355,7 @@ class TMMaxwellOperator(MaxwellOperator): def get_eh_subset(self): return ( (False, False, True) # only ez - + - (True, True, False) # hx and hy + + (True, True, False) # hx and hy ) @@ -400,8 +370,7 @@ class TEMaxwellOperator(MaxwellOperator): def get_eh_subset(self): return ( (True, True, False) # ex and ey - + - (False, False, True) # only hz + + (False, False, True) # only hz ) @@ -416,8 +385,7 @@ class TE1DMaxwellOperator(MaxwellOperator): def get_eh_subset(self): return ( (True, True, False) - + - (False, False, True) + + (False, False, True) ) @@ -432,11 +400,11 @@ class SourceFree1DMaxwellOperator(MaxwellOperator): def get_eh_subset(self): return ( (False, True, False) - + - (False, False, True) + + (False, False, True) ) -def get_rectangular_cavity_mode(E_0, mode_indices): + +def get_rectangular_cavity_mode(E_0, mode_indices): # noqa: N803 """A rectangular TM cavity mode for a rectangle / cube with one corner at the origin and the other at (1,1[,1]).""" dims = len(mode_indices) @@ -452,7 +420,6 @@ def get_rectangular_cavity_mode(E_0, mode_indices): omega = numpy.sqrt(sum(f**2 for f in factors)) - nodes = sym.nodes(dims) x = nodes[0] y = nodes[1] @@ -471,13 +438,13 @@ def get_rectangular_cavity_mode(E_0, mode_indices): tfac = sym.ScalarVariable("t") * omega result = sym.join_fields( - 0, - 0, - sym.sin(kx * x) * sym.sin(ky * y) * sym.cos(tfac), # ez - -ky * sym.sin(kx * x) * sym.cos(ky * y) * sym.sin(tfac) / omega, # hx - kx * sym.cos(kx * x) * sym.sin(ky * y) * sym.sin(tfac) / omega, # hy - 0, - ) + 0, + 0, + sym.sin(kx * x) * sym.sin(ky * y) * sym.cos(tfac), # ez + -ky * sym.sin(kx * x) * sym.cos(ky * y) * sym.sin(tfac) / omega, # hx + kx * sym.cos(kx * x) * sym.sin(ky * y) * sym.sin(tfac) / omega, # hy + 0, + ) else: tdep = sym.exp(-1j * omega * sym.ScalarVariable("t")) diff --git a/grudge/tools.py b/grudge/tools.py index 8a0f5e81e8d0f0a0845163dd0adb40a6ee2560f3..8c5f57b6f0704bc3a5d872cc77ddac5fe366fcf8 100644 --- a/grudge/tools.py +++ b/grudge/tools.py @@ -39,6 +39,8 @@ def is_zero(x): return False +# {{{ SubsettableCrossProduct + class SubsettableCrossProduct: """A cross product that can operate on an arbitrary subsets of its two operands and return an arbitrary subset of its result. @@ -101,6 +103,8 @@ class SubsettableCrossProduct: cross = SubsettableCrossProduct() +# }}} + def count_subset(subset): from pytools import len_iterable diff --git a/setup.cfg b/setup.cfg index b785a55c0d02980068735db0edd18ae81fa70f4a..5a9d9143c76ed52ca9dc4b326cf18d6b1f7f3de4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -3,7 +3,6 @@ ignore = E126,E127,E128,E123,E226,E241,E242,E265,W503,E402 max-line-length=85 exclude= grudge/models/gas_dynamics, - grudge/models/em.py, grudge/models/burgers.py, grudge/models/pml.py, grudge/models/diffusion.py, diff --git a/test/test_grudge.py b/test/test_grudge.py index 288922e07529aa70bffbc1f812efe582e146ff60..e373e1335a86a9f384b19a54bb2a8842130decc8 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -335,7 +335,7 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type @pytest.mark.parametrize("order", [3, 4, 5]) def test_convergence_maxwell(ctx_factory, order): - """Test whether 3D maxwells actually converges""" + """Test whether 3D Maxwell's actually converges""" cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx)