dolfin team mailing list archive

[Evan Lezar <evanlezar@xxxxxxxxx>]

Hi

I am having some trouble implementing driven full-wave problems in dolfin -
most of this from the fact that the field that I need to represent must be a
complex phasor.  I have put together a short document on what it is that I
want to model as well as a small example program that runs, but does not
give the expected results.  It would help immensely if someone could have a
look and tell me if there are any flaws in my reasoning in any of the
stages.

Once this is done, my aim is to improve the handling of complex valued
matrices and functions in dolfin, but it would be useful to have a test case
or two to compare the results to.

Thanks in advance
Evan

Attachment: 3D_waveguide_formulation.pdf
Description: Adobe PDF document

__author__ = "Evan Lezar"
__date__ = "05 Mar 2009"
__copyright__ = "Copyright (C) 2009 Evan Lezar"
__license__  = "GNU LGPL Version 2.1"

"""
    A simple waveguide transmission problem
"""

from dolfin import *
import numpy as N

a = 22.86e-3
b = 10.16e-3
l = 40e-3

k_z_10 = 157.0872
k_o_squared = 43562.7119

def e_xy(x,y):
	return N.sin(N.pi*x/a)

class BoxedDirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        if on_boundary:
            return x[1] > DOLFIN_EPS and abs(x[1] - l) < DOLFIN_EPS
        return False

class Port1(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < DOLFIN_EPS and on_boundary

class Port2(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[1] - l) < DOLFIN_EPS and on_boundary

class U_inc_real(Function):
    def eval(self, values, x):
        values[0] = 0.0
        values[1] = -2*k_z_10*e_xy(x[0],x[1])*N.sin(k_z_10*x[2])
        values[2] = 0
            
class U_inc_imag(Function):
    def eval(self, values, x):
        values[0] = 0.0
        values[1] = -2*k_z_10*e_xy(x[0],x[1])*N.cos(k_z_10*x[2])
        values[2] = 0
        
dolfin_set("linear algebra backend", "uBLAS")
       
mesh = Box(0, 0, 0, a, l, b, 4, 16, 2)
#mesh.refine()

pec = BoxedDirichletBoundary()
port1 = Port1()
port2 = Port2()

sub_domains = MeshFunction("uint", mesh, mesh.topology().dim() - 1)
port1.mark(sub_domains, 1)
port2.mark(sub_domains, 2)

V = FunctionSpace(mesh, "Nedelec", 2)
print "DOFS: ",
print V.dim()

v = TestFunction(V)
w = TrialFunction(V)

DG = FunctionSpace(mesh, "DG", 0)

print 'functions'
#    k_o_squared = Function(DG, '142.33'')
#    k_z = Function(DG, '12.337')

fr = U_inc_real(V)
fi = U_inc_imag(V)
zero = Function(V, '0.0')
print 'done'

volume = 0.5*dot(curl(v), curl(w))*dx - 0.5*k_o_squared*dot(v, w)*dx

n = FacetNormal(mesh)
integrand_real = k_z_10*0.5*dot(cross(n,v), cross(n,w))
surf = integrand_real*ds(1) + integrand_real*ds(2) 

source_real = dot(v, fr)*ds(1)
source_imag = dot(v, fi)*ds(1) 

print "Assembling Matrices"
print "K"
K = assemble(volume)
print "B"
B = assemble(surf, exterior_facet_domains=sub_domains)

print "b"
b_real = assemble(source_real, exterior_facet_domains=sub_domains)
b_imag = assemble(source_imag, exterior_facet_domains=sub_domains)
print "Done"

bc = DirichletBC(V, zero, pec)
bc.apply(K, b_real)
bc.apply(B, b_imag)

print "Done"

A = K.array() + 1j*B.array()
b_source = b_real.array() + 1j*b_imag.array()

print "Solving"
timer = Timer("Solve system")
timer.start()
x = N.linalg.solve(A,b_source)
timer.stop()

print x

field_real = Function(V)
field_real.vector().set(x.real)

field_imag = Function(V)
field_imag.vector().set(x.imag)

e_10 = Function(V, ('0.0', 'sin(pi*x[0]/a)', '0.0'), {'pi':N.pi, 'a':a})

rr = dot(field_real, e_10)
ri = dot(field_imag, e_10)

Rr = assemble(rr*ds(1), exterior_facet_domains=sub_domains)
Ri = assemble(ri*ds(1), exterior_facet_domains=sub_domains)


print "R",
z = 0.0

R = (2*N.exp(-1j*k_z_10*z)/(a*b))*(N.complex(Rr, Ri)) - N.exp(-2j*k_z_10*z)  
print R
print abs(R)

Tr = assemble(rr*ds(2), exterior_facet_domains=sub_domains)
Ti = assemble(ri*ds(2), exterior_facet_domains=sub_domains)

print "T",
z = l
T = 2*N.exp(1j*k_z_10*z)/(a*b)*(N.complex(Tr, Ti))
print T
print abs(T)

summary()