### setting up the objective function

50
views
1
4 weeks ago by
Hello,

I'm simulating the 1D diffusion of a gas. I have experimental data that shows the concentration of the gas near the exit of the system for 1000 time steps. I'm trying to setup a method to minimize the difference between the solution of the PDE at that point and the experimental data point at each time step by changing the diffusion coefficient using dolfin-adjoint. With that being said, I'm not sure if I've set my objective function correctly. I'm also not sure if the solver I'm using is the appropriate one. Any tips, advice or recommended reading would be appreciated (I've looked through the dolfin-adjoint examples online).

Fenics-version:   2017.1.0

Thanks,

Reduced code

from fenics import *
import numpy as np
import pandas as pd
import pyipopt

#The experimental data
colA = temp.ix[:,0]
colB = temp.ix[:,1]

#Time parameters and other parameters for the simulation
T=4.1
num_steps = 1000
dt = T / num_steps
k = Constant(dt)
t_l = 0.036
Diff = 0.5
Diff = Constant(Diff)
Area = 1.1e-5

mesh_size = 20
dz = 1/(mesh_size)
mesh = UnitIntervalMesh(mesh_size)
P1 = FiniteElement('P',mesh.ufl_cell(),1)
V = FunctionSpace(mesh,P1)
v = TestFunction(V)
u = Function(V)
u_n = Function(V)
dx = Measure("dx")

#Function to be solved for (simple diffusion)

#Will have vacuum (C=0) at one end and wall at other (dC/dx = 0)
tol = 1E-14
def boundary_R(x, on_boundary):
return on_boundary and near(x[0],1,tol)
bcs = [DirichletBC(V,Constant(0),boundary_R)]

#Define subdomain so only one point is evaluated (the difference
#between the experimental data point and the value of the function
#at that point)
point_function = MeshFunction("size_t", mesh,0)
point_subdomain = CompiledSubDomain("near(x[0], 0.95)")
point_subdomain.mark(point_function, 1)
dP = Measure('vertex',subdomain_data=point_function)

#The objective function I will try to solve. More will be added with
#each time step. (Unsure if this is the way to set it up or if I'm
#solving what I think I'm solving). colB[0] is a data point(float) at
#the first time step. n will be used later to specify what time step
#it is at
jfunc = assemble(inner(u - colB[0],u - colB[0])*dP(0))

#Perform each iteration of the solution
t = 0
for n in range(num_steps):
t += dt

if t > dt:
pass
else:
str_num = 20
Np = 243.90243902439028
u_L = Np*11500/(Area*t_l)
u_n.vector()[str_num] = u_L

solve(F == 0,u,bcs,solver_parameters={"newton_solver":{"relative_tolerance":1e-10},"newton_solver":{"maximum_iterations":100}})

print("TIME STEP")
print(t)

jfunc += assemble(inner(u - colB[n],u - colB[n])*dP(0))

u_n.assign(u)

#Method for setting up and solving the objective function
control = Control(Diff)
rf = ReducedFunctional(jfunc, control)
problem = MinimizationProblem(rf)
parameters = {"acceptable_tol": 1.0, "maximum_iterations": 1}
solver = IPOPTSolver(problem, parameters=parameters)
solver.solve()
Community: FEniCS Project