Problem with functional

[newbyJulia]

Hello everybody,

i have the following problem with my functional where the inhomogeneous Neumann condition is coupled in the weak formulation.

\int_{\Omega} |\nabla \phi|^2 dx = - \frac{1}{|\Gamma_4} * I_2 * \int_{\Omega_4} \phi ds

I_2 was computed through the Rosenbrock-Wanner method and it is a numpy.ndarray.
When i am trying to optimize the problem with

J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h), dolfin.grad(v))*dx - (leitfaehigkeit*(1/area_gamma1)*((-sol[3,1]))*v*ds(4)))

control = Control(f)

rf = ReducedFunctional(J, control)

problem = MoolaOptimizationProblem(rf)
f_moola = moola.DolfinPrimalVector(f)
solver2 = moola.NewtonCG(problem, f_moola, options={'gtol': 1e-9,
                                                   'maxiter': 20,
                                                   'display': 3,
                                                   'ncg_hesstol': 0})   

sol2 = solver2.solve()

I get the following error in the line solver2.solve()

AttributeError: 'dolfin.cpp.la.Matrix' object has no attribute 'block_variable'

Can somebody please explain me what is wrong?

Greetings

[dokken]

@newbyJulia

\int_{\Omega) |\nabla \phi|^2 dx = - \frac{1}{|\Gamma_4} * I_2 * \int_{\Omega_4} \phi ds

Please write this with Latex formatting, i.e. \int_\Omega (done by $ encapsulation).

Also, please add a minimal code example reproducing hte error, as this makes it alot easier for others to help you, and makes the problem more useful for others that might face the same problem

[newbyJulia]

Thank you dokken. I will do it.

[newbyJulia]

import numpy as np 
import matplotlib.pylab as plt
import math
from numpy.linalg import solve, norm
from dolfin import *
from dolfin_adjoint import *
from mshr import *
import pdb
import moola

# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True


class ode(object):

    # Konstruktor
    def __init__(self, l_1:int, l_2:int, k:float, r_2:float, c_1:float, u_omega:float):
        self.l_1 = l_1
        self.l_2 = l_2
        self.k = k
        self.r_2 = r_2
        self.c_1 = c_1
        self.u_omega = u_omega    

    def rhs(self,q1,p1,q2,p2,l_12):
        f = [p1,
        (-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1))*q1+(l_12/((self.l_1*self.l_2-l_12**2)*self.c_1))*q2,
        p2,
        ((l_12*self.r_2)/((self.l_1*self.l_2-l_12**2)))*p1-((self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2))*p2]
        return f

    def semiImpRK(self,A,q1_i,p1_i,q2_i,p2_i,i,h):        
     # k1
        l_12 = self.k*np.sqrt(self.l_1*self.l_2)     
        b1 = prob.rhs(q1_i,p1_i,q2_i,p2_i,l_12) 
        k1 = np.linalg.solve(A,b1)

        # k2
        b2 = prob.rhs(q1_i+0.5*h*k1[0],p1_i+0.5*h*k1[1],q2_i+0.5*h*k1[2],p2_i+0.5*h*k1[3],l_12)
        k2 = np.linalg.solve(A,b2)

        # yip1
        yip1 = np.zeros(4)
        yip1[0] = q1_i + h*k2[0]
        yip1[1] = p1_i + h*k2[1]
        yip1[2] = q2_i + h*k2[2]
        yip1[3] = p2_i + h*k2[3]

        print("Ergebnis des " + str(i) + "ten Durchlaufs des Algorithmus\n")
        print(str(yip1))
        return yip1[0], yip1[1], yip1[2], yip1[3]

    def problem(self):

        l_12 = self.k*np.sqrt(self.l_1*self.l_2)
    
        jac = np.array([[0, 1, 0, 0],
                    [-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1), 0, l_12/((self.l_1*self.l_2-l_12**2)*self.c_1), 0],
                    [0, 0, 0, 1],
                   [0, (l_12*self.r_2)/(self.l_1*self.l_2-l_12**2), 0, -(self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2)]])

        t0 = 0.0
        T = 100 # Maximale Schweißzeit 15ms
        N = 1000
        h = T/float(N)

        # Die Jacobimatrix muss nur einmal berechnet werden, deshalb liegt diese Berechnung außerhalb der Schleife
        n = 4
        a = 1.0/(2.0+math.sqrt(2.0))
        I = np.identity(n)
        A = I - a*h*jac

        # Anfangswerte
        z = np.zeros((4,N))
        z[0,0] = self.c_1 * self.u_omega
        z[1,0] = 0.0
        z[2,0] = 0.0
        z[3,0] = 0.0

        t = np.zeros(N+1)
        sol = np.zeros((4,N+1))

        # t = 0
        sol[0,0], sol[0,1], sol[0,2], sol[0,3] = prob.semiImpRK(A,z[0,0],z[1,0],z[2,0],z[3,0],0,h) 
        # t > 0
        for i in range(1,N+1):
            q1_i, p1_i, q2_i, p2_i = sol[0,i-1], sol[1,i-1], sol[2,i-1], sol[3,i-1]
            sol[0,i], sol[1,i], sol[2,i], sol[3,i] = prob.semiImpRK(A,q1_i,p1_i,q2_i,p2_i,i,h)
            t[i] = t[i-1] + h

       # Plot-------------------------------------------------
        fig1 = plt.figure(1)
        fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, sharex=True)

        ax1.plot(t, sol[0,:])
        ax1.set_ylabel('A*s')
        ax1.set_title('Q_1')
        ax1.grid(True)

        ax2.plot(t, sol[1,:])
        ax2.set_ylabel('A')
        ax2.set_title('I_1')
        ax2.grid(True)

        ax3.plot(t, sol[2,:])
        ax3.set_ylabel('A*s')
        ax3.set_title('Q_2')
        ax3.grid(True)

        ax4.plot(t, sol[3,:])
        ax4.set_ylabel('A')
        ax4.set_xlabel('Time $t$ [s]')
        ax4.set_title('I_2')
        ax4.grid(True)
        plt.show()
        plt.savefig("Lösung_odeSol.png")

        return sol


class pde(ode):

    def solverIntern(self,a,v,dx,ds,leitfaehigkeit,area_gamma1,sol,f,u_h,bc2,i):

        L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*ds(4)
        dolfin.solve(a == L, u_h, bc2)

        return u_h

    def achsenSetzen(self):
        # Das Achsen-Objekt des Diagramms in einer Variablen ablegen:
        ax = plt.gca()
        # Die obere und rechte Achse unsichtbar machen:
        ax.spines['right'].set_color('none')
        ax.spines['top'].set_color('none')
        # Die linke Diagrammachse auf den Bezugspunkt '0' der x-Achse legen:
        ax.spines['left'].set_position(('data',0))
        # Die untere Diagrammachse auf den Bezugspunkt '0' der y-Achse legen:
        ax.spines['bottom'].set_position(('data',0))
        # Ausrichtung der Achsen-Beschriftung festlegen:
        ax.xaxis.set_ticks_position('bottom')
        ax.yaxis.set_ticks_position('left')

        return ax

    def energyfunctional(self,phi,v,ds,leitfaehigkeit,area_gamma1,mesh,u_h,f,sol):


        W = FunctionSpace(mesh, "DG", 0)
        f = interpolate(Expression("x[0]+x[1]", degree=1), W)

        alpha = Constant(1e-6)
        J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h), dolfin.grad(v))*dx - (leitfaehigkeit*(1/area_gamma1)*((-sol[3,1]))*v*ds(4)))
        print("J = " + str(J))
        control = Control(f)

        rf = ReducedFunctional(J, control)

        problem = MoolaOptimizationProblem(rf)
        f_moola = moola.DolfinPrimalVector(f)
        pdb.set_trace()
        solver2 = moola.NewtonCG(problem, f_moola, options={'gtol': 1e-9,
                                                   'maxiter': 20,
                                                   'display': 3,
                                                   'ncg_hesstol': 0})   


        sol2 = solver2.dolfin_adjoint.solve()
       # sol2 = solver2.solve()
        f_opt = sol['control'].data

        plot(f_opt, title="f_opt")
        plt.savefig("f_opt.png")

        # Define the expressions of the analytical solution
        f_analytic = Expression("1/(1+alpha*4*pow(pi, 4))*w", w=w, alpha=alpha, degree=3)
        u_analytic = Expression("1/(2*pow(pi, 2))*f", f=f_analytic, degree=3)
        #We can then compute the errors between numerical and analytical solutions.

        f.assign(f_opt)
        solve(F == 0, u_h, bc2)
        control_error = errornorm(f_analytic, f_opt)
        state_error = errornorm(u_analytic, u_h)
        print("h(min):           %e." % mesh.hmin())
        print("Error in state:   %e." % state_error)
        print("Error in control: %e." % control_error)

    def __init__(self,l_1,l_2,k,r_2,c_1,u_omega):
        super().__init__(l_1,l_2,k,r_2,c_1,u_omega)

    def solver(self,f,degree,sol):

        fig3 = plt.figure(3)
        # Gebiet und Mesh erstellen------------------------------------------------
        domain_vertices = [Point(2.25, 1.75),
                          Point(1.0, 3.0),
                          Point(-1.0, 1.0),
                          Point(4.0, -4.0),
                          Point(5.0, -3.0),
                          Point(2.0, 0.0),
                          Point(0.5, 0.0)]                


        domain = Polygon(domain_vertices)
        mesh = generate_mesh(domain,20)

        # Plot und Speicherung des nicht konvexen Gebiets---------------------------------
        plt.plot([0.5, 2.25, 1.0, -1.0, 4.0, 5.0, 2.0, 0.5], [0.0, 1.75, 3.0, 1.0, -4.0, -3.0, 0.0, 0.0])
        
        # Das Achsen-Objekt des Diagramms in einer Variablen ablegen:
        ax = prob_pde.achsenSetzen()
        plt.savefig("Domain.png")

        # Plot und Speicherung der Triangulierung------------------------------------------
        plt.figure(4)
        p = plot(mesh, title = 'Finite Element Mesh')
        ax = prob_pde.achsenSetzen()
        plt.savefig("Mesh.png")

        # Randbedingungen-----------------------------------------------------------
        class Left(SubDomain):
            def inside(self, x, on_boundary):
                return on_boundary and (between(x[1], (-0.999, 0.999)) and between(x[0], (3.999, -3.999)))

        class Right(SubDomain):
            def inside(self, x, on_boundary):
                return on_boundary and ((between(x[1], (2.25, 1.75)) and  between(x[0], (1.001, 2.999))) and (between(x[1], (0.5, 0.0)) and between(x[0], (2.249, 1.749))) and (between(x[1], (1.999, 0.0)) and between(x[0], (0.501, 0.0))) and (between(x[1], (4.999, -2.999)) and between(x[0], (2.0, 0.0))))


        class Bottom(SubDomain):
            def inside(self, x, on_boundary):
                return on_boundary and (between(x[1], (5.0, -3.0)) and between(x[0], (4.0, -4.0)))

        class Top(SubDomain):
            def inside(self, x, on_boundary):
                return on_boundary and (between(x[1], (1.0, 3.0)) and between(x[0], (-1.0, 1.0)))

        # Initialize sub-domain instances
        left = Left()
        top = Top()
        right = Right()
        bottom = Bottom()

        # Initialize mesh function for boundary domains
        boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

        boundaries.set_all(4)
        left.mark(boundaries, 1)
        top.mark(boundaries, 2)
        right.mark(boundaries, 3)
        bottom.mark(boundaries, 4)


        # Defining the function space for electric field
        V_phi = FunctionSpace(mesh, "CG", degree)
        #V_phi = VectorFunctionSpace(mesh, "Lagrange", degree)
       # pdb.set_trace()

        # Define Dirichlet boundary condition at top boundary
        bc2 = DirichletBC(V_phi, 0.0, boundaries, 2)
        #Take a look at the solution:
        dx = Measure('dx', mesh)
        ds=Measure('ds', domain=mesh, subdomain_data=boundaries)
        area_gamma1 = dolfin.assemble(1*ds(4))
        leitfaehigkeit = 1.43e-7

        # Permeability (should be 6.3e-3)
        mu = 1e-5

        # Defining the trial function
        phi = TrialFunction(V_phi)

        # Defining the test function
        v = TestFunction(V_phi)

        u_h1 = Function(V_phi, name = "Displacement")

        a = dot(grad(phi), grad(v))*dx

        for i in range(0,10):
            u_h = prob_pde.solverIntern(a,v,dx,ds,leitfaehigkeit,area_gamma1,sol,f,u_h1,bc2,i)           
            p = plot(u_h,title='Finite element solution at t_' + str(i))
            ax = prob_pde.achsenSetzen()
            plot(mesh)
            cbar = plt.colorbar(p)
            cbar.set_label('elektrisches Potential [V]', rotation=270)
            #--------PLOT Gradienten
            plot(-grad(u_h))
            #------------------------
            #plt.savefig('FiniteESol' + str(i)+ '.png')
            plt.savefig("FiniteESol_gmres" + str(i)+ '.png')
            plt.clf()

            # Find the electric field by taking the negative of the
            # gradient of the electric potential. Project this onto
            # the function space for evaluation. (L_2 projection?)
            elec_field = dolfin.project(-grad(u_h))
            plt.figure(5)
            p1 = plot(elec_field)
            cbar = plt.colorbar(p1)
            cbar.set_label('elektrisches Potential [V]', rotation=270)
            plt.show()
            plt.savefig('solution electric_field' + str(i) + '.png')
            plt.clf()

        prob_pde.energyfunctional(phi,v,ds,leitfaehigkeit,area_gamma1,mesh,u_h,f,sol)

    def problem2(self, sol):

        f = Constant(0.0)
        u = prob_pde.solver(f,1,sol)

        print("Test ENDE")

if __name__ == "__main__":

    # Initialisierung einer Instanz der Klasse
    prob = ode(l_1=5.0,
               l_2=1000.0,
               k=0.5,
               r_2=20.0, 
               c_1=8.0,
               u_omega=160)

    prob_pde = pde(l_1=5.0,
               l_2=1000.0,
               k=0.5,
               r_2=20.0, 
               c_1=8.0,
               u_omega=160)
   
    sol = prob.problem()
    prob_pde.problem2(sol)

Maybe it isn’t that bad if you can see everything.

[dokken]

Please reduce your code by removing everything not needed to reproduce the error. This includes plotting etc.

Note that you should try to remove a variable and see if the error is still there, and keep removing variables and reducing complexity of the code until it is a smaller example, such as:

Read before posting: How do I get my question answered?

The following is a list of well constructed questions adhering to the guidelines above: Getting cell dofs for vector valued function space in dolfinx Problem with 3D nedelec element What is the “default” linear solver/preconditioner? Setting Dirichlet boundary condition/constraint on 3D edges

[newbyJulia]

Hello dokken,

here is a minimal working example. I just replaced the Array z with some random numbers.

import numpy as np 
import matplotlib.pylab as plt
import math
from numpy.linalg import solve, norm
from dolfin import *
from dolfin_adjoint import *
from mshr import *
import pdb
import moola


parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

z = np.random.rand(4,1000)

domain_vertices = [Point(2.25, 1.75),
                    Point(1.0, 3.0),
                    Point(-1.0, 1.0),
                    Point(4.0, -4.0),
                    Point(5.0, -3.0),
                    Point(2.0, 0.0),
                    Point(0.5, 0.0)]

domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,20)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (between(x[1], (-0.999, 0.999)) and between(x[0], (3.999, -3.999)))

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((between(x[1], (2.25, 1.75)) and  between(x[0], (1.001, 2.999))) and (between(x[1], (0.5, 0.0)) and between(x[0], (2.249, 1.749))) and (between(x[1], (1.999, 0.0)) and between(x[0], (0.501, 0.0))) and (between(x[1], (4.999, -2.999)) and between(x[0], (2.0, 0.0))))

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (between(x[1], (5.0, -3.0)) and between(x[0], (4.0, -4.0)))

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (between(x[1], (1.0, 3.0)) and between(x[0], (-1.0, 1.0)))

left = Left()
top = Top()
right = Right()
bottom = Bottom()

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

boundaries.set_all(4)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

degree = 2
V_phi = FunctionSpace(mesh, "CG", degree)
bc2 = DirichletBC(V_phi, 0.0, boundaries, 2)
dx = Measure('dx', mesh)
ds=Measure('ds', domain=mesh, subdomain_data=boundaries)
area_gamma1 = dolfin.assemble(1*ds(4))
leitfaehigkeit = 1.43e-7

phi = TrialFunction(V_phi)
v = TestFunction(V_phi)

u_h1 = Function(V_phi, name = "Displacement")

f = Constant(0)
a = dot(grad(phi), grad(v))*dx
for i in range(1,10):
    L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-z[3,i]))*v*ds(4)
    dolfin.solve(a == L, u_h1, bc2)

x = dolfin.SpatialCoordinate(mesh)
W = FunctionSpace(mesh, "DG", 0)
f = interpolate(Expression("x[0]+x[1]", degree=1), W)
J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h1), dolfin.grad(v))*dx) - float(leitfaehigkeit*(1/area_gamma1)*(-z[3,1]))*assemble(v*ds(4))
control = Control(f)

rf = ReducedFunctional(J, control)

problem = MoolaOptimizationProblem(rf)
f_moola = moola.DolfinPrimalVector(f)
solver2 = moola.NewtonCG(problem, f_moola, options={'gtol': 1e-9,
                                                   'maxiter': 20,
                                                   'display': 3,
                                                   'ncg_hesstol': 0})   

sol2 = solver2.solve()
f_opt = sol['control'].data

I am getting the following error:

 func_value = self.scale * self.functional.block_variable.checkpoint
AttributeError: 'dolfin.cpp.la.Vector' object has no attribute 'block_variable'

Greetings

[dokken]

So I do not understand why f is your control, as f is not part of your functional:

@newbyJulia

f = interpolate(Expression("x[0]+x[1]", degree=1), W)
J = assemble(0.5 * dolfin.dot(dolfin.grad(u_h1), dolfin.grad(v))*dx) - float(leitfaehigkeit*(1/area_gamma1)*(-z[3,1]))*assemble(v*ds(4))
control = Control(f)

rf = ReducedFunctional(J, control)

You are also mixing dolfin commands and dolfin_adjoint commands, i.e.

@newbyJulia

for i in range(1,10):
    L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-z[3,i]))*v*ds(4)
    dolfin.solve(a == L, u_h1, bc2)

this solve will not be overloaded, as you are explicitly using the dolfin.solve function.

[newbyJulia]

Hello Mr. Dokken,

thank you for all the good comments about my problem.
I’ve attached a new MWE to ask you to look over it again please.

Of course you were right about the f and the solve. Now i am doing the following:

import numpy as np 
import matplotlib.pylab as plt
import math
from numpy.linalg import solve, norm
from dolfin import *
from dolfin_adjoint import *
from mshr import *
import pdb
import moola
from ufl import nabla_div
import dijitso
import pkgconfig
import scipy.io
import sympy as sy

# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True


class ode(object):

    # Konstruktor
    def __init__(self, l_1:int, l_2:int, k:float, r_2:float, c_1:float, u_omega:float):
        self.l_1 = l_1
        self.l_2 = l_2
        self.k = k
        self.r_2 = r_2
        self.c_1 = c_1
        self.u_omega = u_omega    

    def rhs(self,q1,p1,q2,p2,l_12):
        f = [p1,
        (-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1))*q1+(l_12/((self.l_1*self.l_2-l_12**2)*self.c_1))*q2,
        p2,
        ((l_12*self.r_2)/((self.l_1*self.l_2-l_12**2)))*p1-((self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2))*p2]
        return f

    def semiImpRK(self,A,q1_i,p1_i,q2_i,p2_i,i,h):        
     # k1
        l_12 = self.k*np.sqrt(self.l_1*self.l_2)     
        b1 = prob.rhs(q1_i,p1_i,q2_i,p2_i,l_12) 
        k1 = np.linalg.solve(A,b1)

        # k2
        b2 = prob.rhs(q1_i+0.5*h*k1[0],p1_i+0.5*h*k1[1],q2_i+0.5*h*k1[2],p2_i+0.5*h*k1[3],l_12)
        k2 = np.linalg.solve(A,b2)

        # yip1
        yip1 = np.zeros(4)
        yip1[0] = q1_i + h*k2[0]
        yip1[1] = p1_i + h*k2[1]
        yip1[2] = q2_i + h*k2[2]
        yip1[3] = p2_i + h*k2[3]

        return yip1[0], yip1[1], yip1[2], yip1[3]

    def problem(self):

        l_12 = self.k*np.sqrt(self.l_1*self.l_2)
    
        jac = np.array([[0, 1, 0, 0],
                    [-self.l_2/((self.l_1*self.l_2-l_12**2)*self.c_1), 0, l_12/((self.l_1*self.l_2-l_12**2)*self.c_1), 0],
                    [0, 0, 0, 1],
                   [0, (l_12*self.r_2)/(self.l_1*self.l_2-l_12**2), 0, -(self.l_1*self.r_2)/(self.l_1*self.l_2-l_12**2)]])

        t0 = 0.0
        T = 100 # Maximale Schweißzeit 15ms
        N = 1000
        h = T/float(N)

        # Die Jacobimatrix muss nur einmal berechnet werden, deshalb liegt diese Berechnung außerhalb der Schleife
        n = 4
        a = 1.0/(2.0+math.sqrt(2.0))
        I = np.identity(n)
        A = I - a*h*jac

        # Anfangswerte
        z = np.zeros((4,N))
        z[0,0] = self.c_1 * self.u_omega
        z[1,0] = 0.0
        z[2,0] = 0.0
        z[3,0] = 0.0

        t = np.zeros(N+1)
        sol = np.zeros((4,N+1))

        # t = 0
        sol[0,0], sol[0,1], sol[0,2], sol[0,3] = prob.semiImpRK(A,z[0,0],z[1,0],z[2,0],z[3,0],0,h) 
        # t > 0
        for i in range(1,N+1):
            q1_i, p1_i, q2_i, p2_i = sol[0,i-1], sol[1,i-1], sol[2,i-1], sol[3,i-1]
            sol[0,i], sol[1,i], sol[2,i], sol[3,i] = prob.semiImpRK(A,q1_i,p1_i,q2_i,p2_i,i,h)
            t[i] = t[i-1] + h

        return sol


class pde(ode):

    def __init__(self,l_1,l_2,k,r_2,c_1,u_omega):
        super().__init__(l_1,l_2,k,r_2,c_1,u_omega)

    def solver(self,f,degree,sol):

        mesh = UnitSquareMesh(10, 10)

        class Left(SubDomain):
            def inside(self, x, on_boundary):
                return near(x[0], 0.0)

        class Right(SubDomain):
            def inside(self, x, on_boundary):
                return near(x[0], 1.0)

        class Bottom(SubDomain):
            def inside(self, x, on_boundary):
                return near(x[1], 0.0)

        class Top(SubDomain):
            def inside(self, x, on_boundary):
                return near(x[1], 1.0)

        # Initialize sub-domain instances
        left = Left()
        top = Top()
        right = Right()
        bottom = Bottom()

        # Initialize mesh function for boundary domains
        boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

        boundary_subdomains.set_all(4)
        left.mark(boundary_subdomains, 1)
        top.mark(boundary_subdomains, 2)
        right.mark(boundary_subdomains, 3)
        bottom.mark(boundary_subdomains, 4)

        dss=Measure('ds', domain=mesh, subdomain_data=boundary_subdomains)
        area_gamma1 = dolfin.assemble(1*dss(4))
        # specific electric conductivity of carbon steel at 20 C
        leitfaehigkeit = 1.43e-7

        V_phi = FunctionSpace(mesh, "CG", degree)

        # Dirichlet boundary condition for the equation of potential
        c = Constant(0.0)
        bc2 = DirichletBC(V_phi, c, boundary_subdomains, 2)
        dx = Measure('dx', mesh)

        # Test and trial functions for equation of potential
        phi = TrialFunction(V_phi)
        v = TestFunction(V_phi)

        f = Constant(0.0)

        # time parameters
        T = 15            # final time should be 15ms
        num_steps = 100     # number of time steps
        dt = T / num_steps # time step size
        t = 0.0

        # bilinearform for the equation of potentialß
        a = dot(grad(phi), grad(v))*dx


        # repleasment functions for the equation of potential in the discrete space
        u = Function(V_phi)
        u_n = Function(V_phi)


        for i in range(num_steps):

            # equation of potential not time dependet-------------------------
            L = f*v*dx + (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)
            # Compute solution
            solve(a==L, u, bc2)

            # Update previous solution
            u_n.assign(u)

            t+=dt

            J = assemble(0.5 * dot(grad(phi), grad(v))*dx - (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4))
            f_con = (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)
            control = Control(f_con)
      

    def problem2(self, sol):

        f = Constant(0.0)
        degree = 1
        u = prob_pde.solver(f,degree,sol)

if __name__ == "__main__":

    # Initialisierung einer Instanz der Klasse
    prob = ode(l_1=5.0,
               l_2=1000.0,
               k=0.5,
               r_2=20.0, 
               c_1=8.0,
               u_omega=160)

    prob_pde = pde(l_1=5.0,
               l_2=1000.0,
               k=0.5,
               r_2=20.0, 
               c_1=8.0,
               u_omega=160)
   
    sol = prob.problem()
    prob_pde.problem2(sol)

I am very sorry that I did not remove the ODE class from the MWE, but I think that this is important for my problem.

The error is

AttributeError: 'Form' object has no attribute 'block_variable'.

I want my control to be the voltage drop over the boundary.
May you tell me if there is any possibility for this?

Thank you!

[dokken]

@newbyJulia

   f_con = (leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)
            control = Control(f_con)

You cannot use a whole ufl.Form as a control variable. You need to choose a variable (or create a variable your functional is dependent on, that is your Control).

[newbyJulia]

Thank you for your reply.

Do you mean something like the coordinates of my boundary?
But what if i do not know how this boundary is going to look like after my process?

[dokken]

As your example is far from minimal, it is quite hard for me to help you.
Please formulate the mathematical problem you would like to solve, with functional and PDEs (in strong form).
as of your problem above, you are saying that you would like the integral:

@newbyJulia

(leitfaehigkeit*(1/area_gamma1)*(-sol[3,i]))*v*dss(4)

to be your control. You cannot have an integral as a control, as it depends on several variables.

[newbyJulia]

Oh wow. Thank you.
I have to think about this, because this means i need something what is twice differentiable.

[dokken]

I do not mean that you should implement it. I would simply like to see that mathematical problem you are trying to solve.

[newbyJulia]

Ok.

-\gamma\Delta\phi = 0 \text{ in } \Omega\\ \phi = 0 \text{ on } \Gamma_1\\ -\gamma \frac{\partial \phi}{\partial n} = 0 \text{ on } \Gamma_2 \text{ and } \Gamma_3\\ -\gamma \frac{\partial \phi}{\partial n} = -\frac{1}{|\Gamma_4|}I_2(t) \text{ on } \Gamma_4

whereby $$\phi$$ is the potential, n is the outer normal and $$I_2(t)$$ is a current.
Defining V as

V:=\{v \in H^1|v=0 \text { on } \Gamma_1\}

and multiplying with the test function v and doing the integration yields for the weak formulation

\int_{\Omega} \gamma \nabla \phi \nabla v dx = \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} v(x) dx.

Now we are applying a little trick by saying v should be $$\phi$$ which yields now

\int_{\Omega} \gamma |\nabla \phi|^2 dx = \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} \phi(x) dx.

This would give an optimal control problem with the energy-functional like

\min_{\phi \in V} J(\phi) := \min_{\phi \in V}(\frac{1}{2}\int_{\Omega} \gamma |\nabla \phi |^2 dx - \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} \phi(x) dx)

under the constrains of the potential equation.

I attached a picture of my area and hope that this is understandable (I really can’t draw).

The marked area is to be melted. I was told to only look at the equation of potential and not at the PDEs for the stress and the heat. In an isotropic material J denotes the density of the electric current, E = - \nabla \phi and \phi is the electric potential. With Ohm’s law we have

J = - \gamma(u)\nabla \phi

where u is the temperature.

For the optimal shape of the bulk i should maybe look at the map to the x-axes.

[dokken]

So with this PDE, what is it that you want to use as a control parameter? The current optimal control problem only states \phi as a Control variable, i.e. a classical energy minimization problem.

[newbyJulia]

Yes, this is true.

Indeed i want to minimize the energy, which depends on

\phi.

I think this would be in my code

control = Control(u)

But i do not want to minimize the given energy from the capacitor and at the same time i want to have an optimal shape of my bulk.

[dokken]

This is not the intention of dolfin-adjoint, i would just solve it with the weak form:

@newbyJulia

and multiplying with the test function v and doing the integration yields for the weak formulation

\int_{\Omega} \gamma \nabla \phi \nabla v dx = \frac{1}{|\Gamma_4|}I_2(t)\int_{\Gamma_4} v(x) dx.

∫Ωγ∇ϕ∇vdx=1|Γ4|I2(t)∫Γ4v(x)dx.

[newbyJulia]

Thank you Mr. Dokken.