I want to use Picard iteration to solve the stationary optimal control NS equations.
However I enconter a problem:once I import dolfin_adjoint repository,the function assign() cannot be recognized.
Here is my code:
import matplotlib.pyplot as plt
from dolfin import *
import numpy as np
from dolfin_adjoint import *
#Load mesh and subdomains
mesh = Mesh("data/dolfin_fine.xml")
sub_domains = MeshFunction("size_t", mesh, "data/dolfin_fine_subdomains.xml")
#coefficient
Re=0.01
max_iteration=100
#Define function spaces
V_h = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q_h = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = V_h * Q_h
W = FunctionSpace(mesh, TH)
g=Function(W.sub(0).collapse(),name="control")
P=VectorFunctionSpace(mesh,"Lagrange",2,2)
#No-slip boundary condition for velocity
#x1 = 0, x1 = 1 and around the dolphin
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), g, sub_domains, 0)
#Inflow boundary condition for velocity
#x0 = 1
inflow = Expression(("-sin(x[1]*pi)", "0.0"), degree=2)
bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1)
#Collect boundary conditions
bcs = [bc0, bc1]
#Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
u_n_minus_1=interpolate(f,P)
a = (inner(grad(u), grad(v))+dot(grad(u)*u_n_minus_1,v)- div(v)*p + q*div(u))*dx
L = inner(f, v)*dx
w = Function(W,name="State")
ds=ds(subdomain_data=sub_domains)
n=0
for n in range(max_iteration):
#Update current iteration
print("#################",n+1)
# Compute solution
solve(a == L, w, bcs)
#Split the mixed solution using deepcopy
#(needed for further computation on coefficient vector)
(u, p) = w.split(True)
#Compute error at vertices
vertex_values_u_n_minus_1=u_n_minus_1.compute_vertex_values(mesh)
vertex_values_u=u.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u - vertex_values_u_n_minus_1))
print(error_max)
if error_max<1e-10:
break
#Update previous solution
u_n_minus_1.assign(u)
n+=1
#Split the mixed solution using a shallow copy
(u, p) = w.split()
alpha=Constant(10)
J = assemble(1./2*inner(grad(u), grad(u))*dx + alpha/2*inner(g, g)*ds(0))
m = Control(g)
Jhat = ReducedFunctional(J, m)
g_opt = minimize(Jhat, options={"disp": True})
plot(g_opt, title="Optimised boundary")
g.assign(g_opt)
A,b=assemble_system(a,L,bcs)
solve(A,w.vector(),b)
plot(w.sub(0),title="Velocity")
plot(w.sub(1),title="Pressure")
#Display plots
plt.show()
Here is my error:
Traceback (most recent call last):
File "/home/rsrg/PycharmProjects/pde_project/ns_iterative_bug_algorithm.py", line 80, in <module>
u_n_minus_1.assign(u)
File "/home/rsrg/.local/lib/python3.6/site-packages/fenics_adjoint/types/function.py", line 70, in assign
ret = super(Function, self).assign(other, *args, **kwargs)
File "/usr/lib/python3/dist-packages/dolfin/function/function.py", line 408, in assign
self._cpp_object._assign(rhs._cpp_object)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue发送
*** using the information listed below, you can ask for help at
***
*** [fenics-support@googlegroups.com](mailto:fenics-support@googlegroups.com)
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to initialize vector of degrees of freedom for function.
*** Reason: Cannot re-initialize a non-empty vector. Consider creating a new function.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------