Using solution from FeniCS inside integral

mathematics
1

I’m solving a PDE with Dirac delta initial condition

\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{L_{loc}}\frac{\partial}{\partial\eta}\bigg((\eta^2-1)\frac{\partial p}{\partial\eta}(L,\eta)\bigg), \quad \eta >1,

with initial condition P(L=0) = \delta(\eta-1). With the help of Jeremy B I have code which solves this (\eta has been replaced with x)

"""
1d Fokker-Planck equation, spatially dependant diffusion process with
Dirac delta initial condition, Dirichlet and Neumann.
"""

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt

L = 2          # final L
num_steps = 10   # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term


# Create mesh and define function space
nx = 100
mesh = IntervalMesh(nx, 1., 3)

V = FunctionSpace(mesh, 'P', 1)

x = SpatialCoordinate(mesh)

# Define boundary condition
u_D = Expression("1./(pow(x[0]-y,2) + pow(eps, 2))", eps=1e-6, pi=pi, y=1., degree=2)
# plot(u_D)

def right(x, on_boundary):
    return near(x[0], 1.)
def left(x, on_boundary):
    return near(x[0], 1.)

bc = DirichletBC(V, Constant(0), right)

# Define initial value
u_n = interpolate(u_D, V)
C = assemble(u_n*dx)
u_n.assign(u_n/C)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

a_form = (u*v + (dl/theta)*(x[0]*x[0]-1)*u.dx(0)*v.dx(0))*dx
L_form = u_n*v*dx

# Time-stepping
u = Function(V)
l = 0
for n in range(num_steps):

    # Update current time
    l += dl

    # Compute solution
    solve(a_form == L_form, u, bc)
    
    vtkfile = File('foquesol.pvd')
    vtkfile << u

    # Plot solution
    plot(u)
    #plot(mesh)

    # Update previous solution
    u_n.assign(u)
    print(assemble(u*dx))
   

# Hold plot
plt.xlabel('$\\eta$')
plt.ylabel('$P(\\eta,L)$')
plt.legend(['$P(\\eta,L_0)$','$P(\\eta,L_1)$'], loc='upper right')
plt.show()

How can I use the solution u to evaluate the integral

\int_{1}^{\infty}\bigg(\frac{2}{1+\eta}\bigg)P(L,\eta)d\eta

would I have to perform this in each step of my loop for length stepping?

2

Yes, if you want to compute the above integral for each L you have to do it inside the loop

assemble(2/(1+x[0])*u*dx)
1 Like