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,
Adam

Reduced code

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

#The experimental data
temp = pd.read_csv("../exp_data.txt",delim_whitespace=True,header=None)
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)
F = ((u - u_n) / k)*v*dx  + Diff*dot(grad(u), grad(v))*dx 

#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)

	#Add to the objective function
	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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »