Functional Update through Frequency

dolfin-adjoint
1

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$

Image2
Image22004×135 56.9 KB

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):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

# Load
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'
2

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
3

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)

Instead of looking ahead in time, like:

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

Captura de Tela 2022-05-20 às 09.22.54

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$

Image2
Image22004×135 56.9 KB

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 *
from fenics_adjoint 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):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

# Load

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?

4

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

5

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

6

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.

7

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
8

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?

9

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 *
#from fenics_adjoint 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):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

# Load

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