# Functional Update through Frequency

Dear,

Is it possible to update the Functional through frequency instead of Time?

See the equation below:

$J\left( \Omega \right) =\int^{w_{max}}_{w_{min}} \int_{\Gamma_{N} } [\rho \omega^{2} \left( u^{R}_{next_{i}}w^{R}_{i}\ -\ u^{R}_{i}w^{R}_{i}\ +\ u^{I}_{next_{i}}w^{I}_{i}\ -\ u^{I}_{i}w^{I}_{i}\ \right) dx\ -\ (\sigma (u^{R}_{i}):\nabla w^{R}_{i}\ -\ \sigma (u^{I}_{i}):\nabla w^{I}_{i})dx\ +\ (f^{R}_{i}.w^{R}_{i}\ +\ f^{I}_{i}.w^{I}_{i})ds]dw$


Here’s the code:

from dolfin import *
import numpy

L, H = 2, 1
Nx, Ny = 40, 20
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

#Function space definition over the mesh
Element1 = VectorElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 1)

# Define the mixed element
V_elem = MixedElement([Element1, Element2])

# Define the mixed function space
Vcomplex = FunctionSpace(mesh, V_elem)

#Material properties
E, rho, nu = 1, 1, 0.3

# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
def sigma(v):
dim = v.geometric_dimension()
return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
CompiledSubDomain("x[0]==l && x[1]>=(h/2 - 0.05) && x[1]<=(h/2 + 0.05) ", l=L, h=H).mark(boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
t = Constant((0.0, 1.0))

#Applying the Boundary conditions
zero=Constant((0.0, 0.0))

def left(x, on_boundary):
return near(x[0], 0.)

DirichletBC(Vcomplex.sub(0), zero, bcrigth), DirichletBC(Vcomplex.sub(1), zero, bcrigth)]
boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]

bc=boundary#+source

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
unR, unI = TrialFunctions(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

freqstep = Constant(10)
omega=2*numpy.pi*freqstep

F=omega**2*rho*(dot(unR,wR) - dot(uR,wR) + dot(unI,wI) - dot(uI,wI))*dx - (inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx + (inner(t, wR) + inner(t, wI))*ds

wmin = 0.0 # t
wmax = 100 # end
JtempR = []
JlistR = [JtempR]
JtempI = []
JlistI = [JtempI]
while (wmin <= wmax):

a, L = lhs(F), rhs(F)
solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.solve()
solnUR, solnUI=solnU.split()

unR.assign(solnUR)
unI.assign(solnUI)

wmin += float(freqstep)

JtempR = assemble(inner(t, unR)*ds)
JlistR.append(JtempR)
JtempI = assemble(inner(t, unI)*ds)
JlistI.append(JtempI)

JR = 0
for i in range(1, len(JlistR)):
JR += 0.5*(JlistR[i-1] + JlistR[i])*float(freqstep)

JI = 0
for i in range(1, len(JlistI)):
JI += 0.5*(JlistI[i-1] + JlistI[i])*float(freqstep)


Error:

unR.assign(solnUR)

AttributeError: 'ListTensor' object has no attribute 'assign'


I don’t understand the question, although the error can be avoided with the following modifications:

#unR, unI = TrialFunctions(Vcomplex)
Un = Function(Vcomplex)


and

    #unR.assign(solnUR)
#unI.assign(solnUI)
Un.assign(solnU)


There is also another error at the end, which can be avoided by

#JtempR = []
#JlistR = [JtempR]
JlistR = []
#JtempI = []
#JlistI = [JtempI]
JlistI = []


although it’s difficult for me to say whether that’s “correct”, as I’m not sure what the intended behavior of the code is.

In any case, it looks like this is a workaround MixedElement implementation of complex arithmetic. This could probably be simplified a lot by instead using the native complex number support in FEniCSx (circumventing the need for a MixedElement and simplifying the variational form), although I don’t work with complex-valued PDEs much myself, and haven’t tried using the new complex number support.

1 Like

Dear @kamensky, thank you very much for your support.

The idea is to carry out the analysis in the frequency domain by varying the frequency itself. That’s why I’m asking you if what I’m “thinking” about doing numerically is right. I’m following the time parsing guidelines as well: Expressing functionals

And that’s still the reason I don’t go straight to FEniCSx, because of the dolfin-adoint (which I think is very elegant and I don’t understand why it seems that there are few users and support here)

$J(u) = \int_0^T\int_{\Omega}\left\langle u(t),u(t)\right\rangle \ \textrm{d}\Omega \ \textrm{d}t$


Look at the frequency advance, such as:

$J\left( \Omega \right) =\int^{w_{max}}_{w_{min}} \int_{\Gamma_{N} } [\rho \omega^{2} \left( u^ {R}_{next_{i}}w^{R}_{i}\ -\ u^{R}_{i}w^{R}_{i}\ +\ u^{I}_{ next_{i}}w^{I}_{i}\ -\ u^{I}_{i}w^{I}_{i}\ \right) dx\ -\ (\sigma (u^{ R}_{i}):\nabla w^{R}_{i}\ -\ \sigma (u^{I}_{i}):\nabla w^{I}_{i})dx\ +\ (f^{R}_{i}.w^{R}_{i}\ +\ f^{I}_{i}.w^{I}_{i})ds]dw$


Assuming wmin is 0.

(I don’t know why the formulas aren’t getting right, I’m putting it in between  so I’m putting the images as well.)

And with that, as per your directions, I changed the code, and it looks like it worked!

See the code:

from dolfin import *
import numpy

L, H = 2, 1
Nx, Ny = 40, 20
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

#Function space definition over the mesh
Element1 = VectorElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 1)

# Define the mixed element
V_elem = MixedElement([Element1, Element2])

# Define the mixed function space
Vcomplex = FunctionSpace(mesh, V_elem)

#Material properties
E, rho, nu = 1, 1, 0.3

# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
def sigma(v):
dim = v.geometric_dimension()
return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
CompiledSubDomain("x[0]==l && x[1]>=(h/2 - 0.05) && x[1]<=(h/2 + 0.05) ", l=L, h=H).mark(boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
t = Constant((0.0, 1.0))

#Applying the Boundary conditions
zero=Constant((0.0, 0.0))

def left(x, on_boundary):
return near(x[0], 0.)

boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]

bc=boundary

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
unR, unI = TrialFunctions(Vcomplex)
Un = Function(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

freqstep = Constant(10)

omega=2*numpy.pi*freqstep

F=omega**2*rho*(dot(unR,wR) - dot(uR,wR) + dot(unI,wI) - dot(uI,wI))*dx - (inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx + (inner(t, wR) + inner(t, wI))*ds

wmin = 0.0 # t
wmax = 100 # end
#JtempR = []
#JlistR = [JtempR]
JlistR = []
#JtempI = []
#JlistI = [JtempI]
JlistI = []
while (wmin <= wmax):

a, L = lhs(F), rhs(F)
solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.solve()

solnUR, solnUI=solnU.split()

Un.assign(solnU)
#unR.assign(solnUR)
#unI.assign(solnUI)

wmin += float(freqstep)

JtempR = assemble(inner(t, solnUR)*ds)
JlistR.append(JtempR)
JtempI = assemble(inner(t, solnUI)*ds)
JlistI.append(JtempI)

JR = 0
for i in range(1, len(JlistR)):
JR += 0.5*(JlistR[i-1] + JlistR[i])*float(freqstep)

JI = 0
for i in range(1, len(JlistI)):
JI += 0.5*(JlistI[i-1] + JlistI[i])*float(freqstep)


From a mathematical and numerical point of view, is this correct?

I actually believe it’s still not working. If you look at JlistR and JlistI they have the same value, it’s not updating.

Dear @kamensky, Do you understand what I want to do, could you help me?

I’m still not sure what the mathematical problem statement or intended behavior of the code are. For instance, the following points are still unclear to me: Is \omega supposed to be the integration variable in the outer integral? (I.e., should dw be d\omega instead?) Does the functional J involve a derivative with respect to \omega/w that is being approximated as a finite difference? Is there maybe some confusion in the equations and/or code over the functional J in the equations vs. the PDE residual F in the code?

A quick note on formatting equations here is that you only need dollar signs to typeset equations with LaTeX syntax; there doesn’t need to be an extra layer of backticks around the dollar signs.

If I understand correctly, this is a forced harmonic analysis, and the intent is to parametrically vary the excitation frequency, omega. In that case, you need to update the value of omega after incrementing wmin, i.e.

wmin += float(freqstep)
omega.assign(wmin)


Two other notes:

1. Your variational formulation appears to be incorrect. You have redundant definitions of uR, uI = TrialFunctions(Vcomplex) and unR, unI = TrialFunctions(Vcomplex). The term
omega**2*rho*(dot(unR,wR) - dot(uR,wR) + dot(unI,wI) - dot(uI,wI))*dx

is equal to 0 (because unR == uR and unI == uI, so the omega-dependence is removed.
2. You should consider moving the definition of problem and solver outside of the loop to reduce memory usage.

By changing the last part of your code to the following, the values in JlistR and JlistI update on each step.

freqstep = 10.0
freq = Constant(freqstep)
omega=2*numpy.pi*freq

F=omega**2*rho*(dot(unR,wR) + dot(unI,wI))*dx - (inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx + (inner(t, wR) + inner(t, wI))*ds

wmin = 0.0 # t
wmax = 10*freqstep # end
#JtempR = []
#JlistR = [JtempR]
JlistR = []
#JtempI = []
#JlistI = [JtempI]
JlistI = []

a, L = lhs(F), rhs(F)
solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
while (wmin <= wmax):

solver.solve()

solnUR, solnUI=solnU.split()

Un.assign(solnU)
#unR.assign(solnUR)
#unI.assign(solnUI)

wmin += float(freqstep)
freq.assign(wmin)

JtempR = assemble(inner(t, solnUR)*ds)
JlistR.append(JtempR)
JtempI = assemble(inner(t, solnUI)*ds)
JlistI.append(JtempI)

1 Like

Dear @conpierce8,

That’s exactly what I’m trying to implement. Thank you!

Well, however, if you insert at the end an analysis to perform a time integral numerically, for example, by the trapezoidal rule:

JR = 0
for i in range(1, len(JlistR)):
JR += 0.5*(JlistR[i-1] + JlistR[i])*float(freqstep)

JI = 0
for i in range(1, len(JlistI)):
JI += 0.5*(JlistI[i-1] + JlistI[i])*float(freqstep)


There is no JR or JI update, do you know why?

Can you clarify what you mean by “no JR or JI update”? When I add some print statements to the code, it appears as though JR and JI are updating. Here’s my code and the output from the code:

from dolfin import *
import numpy

L, H = 2, 1
Nx, Ny = 40, 20
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

#Function space definition over the mesh
Element1 = VectorElement("CG", mesh.ufl_cell(), 1)
Element2 = VectorElement("CG", mesh.ufl_cell(), 1)

# Define the mixed element
V_elem = MixedElement([Element1, Element2])

# Define the mixed function space
Vcomplex = FunctionSpace(mesh, V_elem)

#Material properties
E, rho, nu = 1, 1, 0.3

# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
def sigma(v):
dim = v.geometric_dimension()
return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
CompiledSubDomain("x[0]==l && x[1]>=(h/2 - 0.05) && x[1]<=(h/2 + 0.05) ", l=L, h=H).mark(boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
t = Constant((0.0, 1.0))

#Applying the Boundary conditions
zero=Constant((0.0, 0.0))

def left(x, on_boundary):
return near(x[0], 0.)

boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]

bc=boundary

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
unR, unI = TrialFunctions(Vcomplex)
Un = Function(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

freqstep = 10.0
freq = Constant(freqstep)
omega=2*numpy.pi*freq

F=omega**2*rho*(dot(unR,wR) + dot(unI,wI))*dx - (inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx + (inner(t, wR) + inner(t, wI))*ds

wmin = 0.0 # t
wmax = 10*freqstep # end
#JtempR = []
#JlistR = [JtempR]
JlistR = []
#JtempI = []
#JlistI = [JtempI]
JlistI = []

a, L = lhs(F), rhs(F)
solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
while (wmin <= wmax):

solver.solve()

solnUR, solnUI=solnU.split()

Un.assign(solnU)
#unR.assign(solnUR)
#unI.assign(solnUI)

wmin += float(freqstep)
freq.assign(wmin)

JtempR = assemble(inner(t, solnUR)*ds)
JlistR.append(JtempR)
JtempI = assemble(inner(t, solnUI)*ds)
JlistI.append(JtempI)

JR = 0
for i in range(1, len(JlistR)):
JR += 0.5*(JlistR[i-1] + JlistR[i])*float(freqstep)
print("i =", i, "JR =", JR)

JI = 0
for i in range(1, len(JlistI)):
JI += 0.5*(JlistI[i-1] + JlistI[i])*float(freqstep)
print("i =", i, "JI =", JI)

i = 1 JR = -0.977993157707851
i = 2 JR = -1.608793567237851
i = 3 JR = -1.8042159773626132
i = 4 JR = -1.8868042562189538
i = 5 JR = -1.933995643984482
i = 6 JR = -1.9647563711541776
i = 7 JR = -1.9864581207037124
i = 8 JR = -2.00261173706925
i = 9 JR = -2.015113292078411
i = 10 JR = -2.0250800708636785
i = 1 JI = -0.9779931577081202
i = 2 JI = -1.6087935672382547
i = 3 JI = -1.8042159773630169
i = 4 JI = -1.8868042562193574
i = 5 JI = -1.9339956439848858
i = 6 JI = -1.9647563711545812
i = 7 JI = -1.986458120704116
i = 8 JI = -2.0026117370696537
i = 9 JI = -2.015113292078815
i = 10 JI = -2.025080070864082

1 Like