Optimal Control Poisson Equation on Subdomain

dolfin-adjoint
1

Hey guys,

I’m a little bit frustrated. I have following code which solves an optimal control problem with the poisson equation:

from dolfin import *
from dolfin_adjoint import *
import moola
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

set_log_level(LogLevel.ERROR)

# Speicher für Ergebnisse
results = []

# Verschiedene Gittergrößen
for n in [32, 64, 128, 256]:
    print(f"Löse für n = {n} ...")
    
    mesh = UnitSquareMesh(n, n)

    # Lokale Verfeinerung
   # cf = MeshFunction("bool", mesh, mesh.geometric_dimension())
   # subdomain = CompiledSubDomain("std::abs(x[0]-0.5) < 0.25 && std::abs(x[1]-0.5) < 0.25")
   # subdomain.mark(cf, True)
   # mesh = refine(mesh, cf)

    # Funktionräume
    V = FunctionSpace(mesh, "CG", 1)
    W = FunctionSpace(mesh, "DG", 0)

    f = Function(W, name="Control")
    u = Function(V, name="State")
    v = TestFunction(V)

    # Zustandsgleichung
    F = (inner(grad(u), grad(v)) - f * v) * dx
    bc = DirichletBC(V, 0.0, "on_boundary")
    solve(F == 0, u, bc)

    # Steuerung & Ziel
    x = SpatialCoordinate(mesh)
    w = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=3)
    d = 1 / (2 * pi ** 2)
    d = Expression("d*w", d=d, w=w, degree=3)

    alpha = Constant(1e-6)
    J = assemble((0.5 * inner(u - d, u - d)) * dx + alpha / 2 * f ** 2 * dx)
    control = Control(f)
    
    # reduziertes Zielfunktional
    rf = ReducedFunctional(J, Control(f))

    # Optimierungsproblem lösen mit CG Verfahren
    problem = MoolaOptimizationProblem(rf)
    f_moola = moola.DolfinPrimalVector(f)
    solver = moola.NewtonCG(problem, f_moola, options={
        'gtol': 1e-9,
        'maxiter': 20,
        'display': 0,
        'ncg_hesstol': 0
    })

    sol = solver.solve()
    f_opt = sol["control"].data
   
    # Definition der analytischen Lösung um die Fehler zu berechnen
    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)


    f_analytic_v = interpolate(f_analytic, W)
    f.assign(f_opt)

    # Zustand mit optimierter Steuerung berechnen
    solve(F == 0, u, bc)

    # Fehlerberechnung
    control_error = errornorm(f_analytic, f_opt)
    state_error = errornorm(u_analytic, u)
    h_min = mesh.hmin()

    # Ergebnis speichern
    results.append((n, h_min, state_error, control_error))

# Datenframe erstellen
df = pd.DataFrame(results, columns=["n", "h_min", "state_error", "control_error"])

# Konvergenzraten berechnen
def convergence_rate(errors, hs):
    return [None] + [np.log(errors[i-1]/errors[i]) / np.log(hs[i-1]/hs[i]) for i in range(1, len(errors))]

df["convergence_rate_state"] = convergence_rate(df["state_error"], df["h_min"])
df["convergence_rate_control"] = convergence_rate(df["control_error"], df["h_min"])

# Tabelle anzeigen
print("Konvergenztabelle:")
print(df.to_string(index=False, float_format="%.2e"))

# Plot
plt.figure(figsize=(8, 6))
plt.loglog(df["h_min"], df["state_error"], "o-", label="State-Fehler")
plt.loglog(df["h_min"], df["control_error"], "s--", label="Control-Fehler")
plt.loglog(df["h_min"], df["state_error"].iloc[0]*(df["h_min"]/df["h_min"].iloc[0])**2, "k:", label="O(h²)")
plt.loglog(df["h_min"], df["control_error"].iloc[0]*(df["h_min"]/df["h_min"].iloc[0])**1, "k-.", label="O(h)")
plt.xlabel(r"$h_{\min}$")
plt.ylabel("Fehler")
plt.title("Konvergenzverhalten (Log-Log)")
plt.grid(True, which="both", ls=":")
plt.legend()
plt.tight_layout()
plt.show()



plot(f, title="Optimale Steuerung (f)")

plt.show()

This works fine. But now I am trying to constrain my control on a subdomain, so it only works there and is == 0 outside of the subdomain. I tried to update my code with this class:

from dolfin import *
from dolfin_adjoint import *
import moola
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

set_log_level(LogLevel.ERROR)

# Speicher für Ergebnisse
results = []

# Verschiedene Gittergrößen
for n in [32, 64, 128, 256]:
    print(f"Löse für n = {n} ...")
    
    mesh = UnitSquareMesh(n, n)
    class ControlRegion(SubDomain):
        def inside(self, x, on_boundary):
            return (0.25 < x[0] < 0.75) and (0.25 < x[1] < 0.75)

    subdomain_markers = MeshFunction("size_t", mesh, mesh.topology().dim())
    subdomain_markers.set_all(0)
    omega = ControlRegion()
    omega.mark(subdomain_markers, 1)
    dx_sub = Measure("dx", domain=mesh, subdomain_data=subdomain_markers)
    plot(subdomain_markers)
    plt.title("Control Region (omega)")
    plt.show()

    # Lokale Verfeinerung
   # cf = MeshFunction("bool", mesh, mesh.geometric_dimension())
   # subdomain = CompiledSubDomain("std::abs(x[0]-0.5) < 0.25 && std::abs(x[1]-0.5) < 0.25")
   # subdomain.mark(cf, True)
   # mesh = refine(mesh, cf)

    # Funktionräume
    V = FunctionSpace(mesh, "CG", 1)
    W = FunctionSpace(mesh, "DG", 0)

    f = Function(W, name="Control")
    u = Function(V, name="State")
    v = TestFunction(V)

    # Zustandsgleichung
    F = (inner(grad(u), grad(v)) - f * v) * dx
    bc = DirichletBC(V, 0.0, "on_boundary")
    solve(F == 0, u, bc)

    # Steuerung & Ziel
    x = SpatialCoordinate(mesh)
    w = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=3)
    d = 1 / (2 * pi ** 2)
    d = Expression("d*w", d=d, w=w, degree=3)

    alpha = Constant(1e-6)
    J = 0.5 * inner(u - d, u - d) * dx + alpha / 2 * f**2 * dx_sub(1)

    control = Control(f)
    
    # reduziertes Zielfunktional
    rf = ReducedFunctional(J, Control(f))

    # Optimierungsproblem lösen mit CG Verfahren
    problem = MoolaOptimizationProblem(rf)
    f_moola = moola.DolfinPrimalVector(f)
    solver = moola.NewtonCG(problem, f_moola, options={
        'gtol': 1e-9,
        'maxiter': 20,
        'display': 0,
        'ncg_hesstol': 0
    })

    sol = solver.solve()
    f_opt = sol["control"].data
   
    # Definition der analytischen Lösung um die Fehler zu berechnen
    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)


    f_analytic_v = interpolate(f_analytic, W)
    f.assign(f_opt)

    # Zustand mit optimierter Steuerung berechnen
    solve(F == 0, u, bc)

    # Fehlerberechnung
    control_error = errornorm(f_analytic, f_opt)
    state_error = errornorm(u_analytic, u)
    h_min = mesh.hmin()

    # Ergebnis speichern
    results.append((n, h_min, state_error, control_error))

# Datenframe erstellen
df = pd.DataFrame(results, columns=["n", "h_min", "state_error", "control_error"])

# Konvergenzraten berechnen
def convergence_rate(errors, hs):
    return [None] + [np.log(errors[i-1]/errors[i]) / np.log(hs[i-1]/hs[i]) for i in range(1, len(errors))]

df["convergence_rate_state"] = convergence_rate(df["state_error"], df["h_min"])
df["convergence_rate_control"] = convergence_rate(df["control_error"], df["h_min"])

# Tabelle anzeigen
print("Konvergenztabelle:")
print(df.to_string(index=False, float_format="%.2e"))

# Plot
plt.figure(figsize=(8, 6))
plt.loglog(df["h_min"], df["state_error"], "o-", label="State-Fehler")
plt.loglog(df["h_min"], df["control_error"], "s--", label="Control-Fehler")
plt.loglog(df["h_min"], df["state_error"].iloc[0]*(df["h_min"]/df["h_min"].iloc[0])**2, "k:", label="O(h²)")
plt.loglog(df["h_min"], df["control_error"].iloc[0]*(df["h_min"]/df["h_min"].iloc[0])**1, "k-.", label="O(h)")
plt.xlabel(r"$h_{\min}$")
plt.ylabel("Fehler")
plt.title("Konvergenzverhalten (Log-Log)")
plt.grid(True, which="both", ls=":")
plt.legend()
plt.tight_layout()
plt.show()



plot(f, title="Optimale Steuerung (f)")

plt.show()

The plot of the subdomain markers shows exactly what I want but I can’t get it to run through the optimization.
Does somebody maybe have a similiar problem and can assisst me?
Thanks in advance