Dear community,
The minimum working code example is based on the ft03_heat.py tutorial, which describes a time dependent poisson equation and is governed by the following equations:
u’=Δu+f in the unit square
u=uD on the boundary (Dirichlet boundary conditions)
u=u0 at t=0
I tried to implement the time-, space and u-dependend source term f by the following code
# Define source term f
class MySourceConcentration(Expression):
def __init__(self, t, u):
self.t = None
self.u = None
def eval(self, value, t, x, u, area_old):
# Calculate the source term f for t>0
if t>0:
if x[0] >= 0.0 and x[0] <= 0.1:
# Calculate average of u for the domain x>=0 and x<=0.1
u_avg = assemble(u*dx) / assemble(Constant(1.0)*dx)
# Calculate the surface area of a particle (equation is made up for demonstrationen purposes without physical correctness)
area = area_0 - dt * sqrt(aera_old) * (1 - u_avg)
# Calculate source term f
value[0] = Constant1 * area * (Constant2 - u_avg)
# set area to area old for next time step
area_old = area
# for the rest of the domain: source term f=0
else:
value[0] = 0
return value[0]
# Set the source term f for t=0 to initial value f_0 (to provide a physical sound example, ...
# ...there has also to be a differentiation regarding the space, but to reduce the complexity of the example, I just set f=f_0 for the whole domain)
else:
value[0]=f_0
return value[0]
def value_shape(self):
return ()
and called the source term f by the following code line
f = MySourceConcentration(t, u_n)
But no matter how I change the code, I get a error message similar to this one:
Traceback (most recent call last):
File "20200416_ft03_heat_f_class.py", line 94, in <module>
f = MySourceConcentration(t, u_n)
File "20200416_ft03_heat_f_class.py", line 7, in __init__
self.t = 0
File "/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py", line 438, in __setattr__
elif name in self._parameters:
File "/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py", line 432, in __getattr__
return self._parameters[name]
File "/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py", line 432, in __getattr__
return self._parameters[name]
File "/usr/local/lib/python3.6/dist-packages/dolfin/function/expression.py", line 432, in __getattr__
return self._parameters[name]
[Previous line repeated 327 more times]
RecursionError: maximum recursion depth exceeded
After reading the Python documentation to classes and checking the Fenics book and not finding a solution, I thought it’s time to post a thread on the fenics discourse group.
On top of that, I´m not sure if
# set area to area old for next time step
area_old = area
of the eval function is appropriated to set area_old for the next time step or in other words, I´m not sure, if the eval function will remember area_old at the beginning of the new time step to calculate area by
area = area_0 - dt * sqrt(aera_old) * (1 - u_avg)
Thanks a lot!
Click to see my full MWCE.
from fenics import *
import numpy as np
# Define source term f
class MySourceConcentration(Expression):
def __init__(self, t, u):
self.t = None
self.u = None
def eval(self, value, t, x, u, area_old):
# Calculate the source term f for t>0
if t>0:
if x[0] >= 0.0 and x[0] <= 0.1:
# Calculate average of u for the domain x>=0 and x<=0.1
u_avg = assemble(u*dx) / assemble(Constant(1.0)*dx)
# Calculate the surface area of a particle (equation is made up for demonstrationen purposes without physical correctness)
area = area_0 - dt * sqrt(aera_old) * (1 - u_avg)
# Calculate source term f
value[0] = Constant1 * area * (Constant2 - u_avg)
# set area to area old for next time step
area_old = area
# for the rest of the domain: source term f=0
else:
value[0] = 0
return value[0]
# Set the source term f for t=0 to initial value f_0 (to provide a physical sound example, ...
# ...there has also to be a differentiation regarding the space, but to reduce the complexity of the example, I just set f=f_0 for the whole domain)
else:
value[0]=f_0
return value[0]
def value_shape(self):
return ()
T = 2.0 # final time
num_steps = 10 # number of time steps
dt = T / num_steps # time step size
alpha = 3 # parameter alpha
beta = 1.2 # parameter beta
Constant1 = 1
Constant2 = 10
C_S = 5
area_0 = 5
area_old = area_0
f_0 = 4
# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'P', 1)
# Create classes for defining parts of the boundaries and the interior
# of the domain
class Obstacle_left(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (0.0, 0.1)))
# Initialize sub-domain instances
obstacle_left = Obstacle_left()
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, space_dim-1)
obstacle_left.mark(domains, 1)
# Define new measures associated with the interior domains
dx = Measure('dx', domain=mesh, subdomain_data=domains)
vtkfile = File("20200416_ft03_heat_f_class/solution.pvd")
# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
degree=2, alpha=alpha, beta=beta, t=0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define initial value
u_n = interpolate(u_D, V)
# Define variational problem
t=0
u = TrialFunction(V)
v = TestFunction(V)
f = MySourceConcentration(t, u_n)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f[0])*v*dx
a, L = lhs(F), rhs(F)
# Time-stepping
u = Function(V)
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(a == L, u, bc)
vtkfile << (u, t)
# Update previous solution
u_n.assign(u)
Python to use