Solving PDE with First and Second Space Derivatives

[Esteban Lopez]

Hello,

I am currently attempting to solve the following PDE, with the following initial and boundary conditions using pyomo.dae.

PDE

Initial Condition

Boundary Condition Left Side

Boundary Condition Right Side

This is a typical one-dimensional advection equation for solving a mass transfer problem in a pipe of length L and constant velocity u, and an axial dispersion term Dax. It includes first and second order space derivatives.

While I am able to set up the optimization model in pyomo and solve it using the ASL interface, with ipopt_sens for the non-linear solver and ma27 for the linear solver.

The model solution seems to work fine for low Dax values. However, for high Dax values, even though the optimal solution is found, I keep getting a wrong solution.

This is a well-studied problem, so I know beforehand how the solution is supposed to look for high Dax values, and that's how I know I'm getting the incorrect solution. This is my code that sets up the model and calls the solvers. Any help will be really appreciated. 

I have been able to solve this problem in Matlab and get the correct solution for all types of Dax values. I am a starter in Pyomo, so I'm not sure if I'm going wrong in setting up the model itself, the discretization or calling the solver (or something else). With the values of L and u from the code values, I start seeing problem at Peclet numbers smaller than 100. Peclet number is  Pe = (L*U)/Dax. 

This is a picture of the solutions at various Peclet number s using Matlab for reference (Matlab uses a finite difference for both space and time, and I use 100 finite elements for each time and space).

And this is a pic from pyomo results.

Thank you,

Esteban

from pyomo.environ import *

from pyomo.dae import *

from pyomo.opt import SolverFactory

import matplotlib.pyplot as plt

# -------------------------------------------------------------------

####Define Parameters and Constants

# -------------------------------------------------------------------

#Inital Concentration

C0_gL= 8 #g/L

C0= C0_gL/1000

#Length of the column cm

L=20

#Time of operation hr

tend=3

#Diamter cm

d=0.94

#Bed Area

Ab=np.pi*d**2/4

#Flow rate cc/min

q=0.2

#Porosity

eb=0.539

#Velocity (u-> interstitial velocity) cm/h

u_superficial=q/Ab*60

u=u_superficial/eb

#Peclet number vL/Dax

Pe=1000.

#Axial Dispersion Coefficient cm^2/h

Dax=u*L/Pe

#Number of discretization points in time t and space x

nfet=30

nfex=10

#Number of collocation points if the collocation points method is used

ncp=3

m=ConcreteModel()

# -------------------------------------------------------------------

####Define PYOMO Model

# -------------------------------------------------------------------

#Define model parameters as pyomo params

m.C0 = Param(initialize=C0)

m.Dax= Param(initialize = Dax)

m.L = Param(initialize = L)

m.eb = Param(initialize = eb)

m.U = Param(initialize=u)

#Define Continous set of independent variables

m.x = ContinuousSet(bounds=(0.,1.))

m.t = ContinuousSet(bounds=(0.,1.))

#Define Variables and other parameters

m.C = Var(m.t, m.x)

#Define derivative variables

m.dCdt = DerivativeVar(m.C, wrt=m.t) # differentiated liquid concentration respect to time

m.dCdx = DerivativeVar(m.C, wrt=m.x)

m.dC2dx2 = DerivativeVar(m.C, wrt=(m.x, m.x))

# -------------------------------------------------------------------

####Define Constraints

# -------------------------------------------------------------------

#Constraints

#Mass Balance Constraint

def MassBalanceLiquid_rule(m, k,l):

    if (k==0. or l==0. or l==1.): return Constraint.Skip

    else: return m.dCdt[k,l]/tend + m.U*m.dCdx[k,l]/m.L -m.Dax*m.dC2dx2[k,l]/m.L/m.L == 0

    

#Write up the constraints

m.MassBalanceLiquid=Constraint(m.t, m.x, rule=MassBalanceLiquid_rule)

# -------------------------------------------------------------------

####Define BC0s and IC

# -------------------------------------------------------------------

#Boundary Condition, at the entrance of the column (x=0) concentration is C0, for all t>0

#Initial Condition, C = 0 for all x

def InitialCondition_rule(m, l):

    if l==0 or l==1.: return Constraint.Skip

    else: return m.C[0.,l]==0

def BoundaryCondition_rule(m, k):

    return m.U*m.C0==m.U*m.C[k,0.]-m.Dax*m.dCdx[k,0.]/m.L

      

def BoundaryCondition_rule2(m, k):

    return m.dCdx[k,1.]/m.L==0

#Write up these conditions as constraints   

m.InitialCondition=Constraint(m.x, rule=InitialCondition_rule)

m.BoundaryCondition=Constraint(m.t, rule=BoundaryCondition_rule)

m.BoundaryCondition2=Constraint(m.t, rule=BoundaryCondition_rule2)

# -------------------------------------------------------------------

#### Discretization

# -------------------------------------------------------------------

#Set up the discretization grid, use 'Backward' for finite_difference

#Use finite difference for x and collocation for t

#time discretization

discretizet=TransformationFactory('dae.collocation')

discretizet.apply_to(m,nfe=nfet, ncp=ncp , wrt=m.t, scheme='LAGRANGE-RADAU')

#x discretization

xscheme='BACKWARD'

discretizex=TransformationFactory('dae.finite_difference')

discretizex.apply_to(m,nfe=nfex,wrt=m.x,scheme=xscheme)

# -------------------------------------------------------------------

# Solver options

# -------------------------------------------------------------------

#### Create the ipopt_sens solver plugin using the ASL interface

solver = 'ipopt_sens'

solver_io = 'nl'

stream_solver = True    # True prints solver output to screen

keepfiles =     False    # True prints intermediate file names (.nl,.sol,...)

opt = SolverFactory(solver)#,solver_io=solver_io)

####

opt.options['mu_init'] = 1e-3

#opt.options['ma57_pivtol'] = 1e-8

opt.options['max_iter'] = 5000

#opt.options['linear_system_scaling'] = 'mc19'

#

# opt.options['linear_solver'] = 'ma97'

# opt.options['linear_solver'] = 'ma57'

opt.options['linear_solver'] = 'ma27'

#m.preprocess()

results = opt.solve(m, tee=True)

#######################Plotting

#Collect Results for Plotting

r=len(m.t)

c=len(m.x)

C_num=np.zeros(r)

t=np.zeros(r)

for i in range(r): t[i]=m.t.card(i+1)

for j in range(r):

    C_num[j]=value(m.C[m.t.card(j+1),m.x.card(-1)])

#Plot

plt.rcParams.update({'font.size': 14})

plt.plot(t*tend,C_num/C0,label='Pe={:}'.format(Pe))

plt.ylabel('C/C0')

plt.xlabel('Time,t (hr)')

plt.ylim((0,1))

plt.legend()

plt.show()

[Esteban Lopez]

I actually found the answer to my problem. Realized I had to define the first order space derivative at x=0, since I was using a backward finite difference scheme.