### [Solved]Trouble met when trying to simulate a model with changing parameters.

157

views

1

```
# coding: utf-8
# In[ ]:
"""
FEniCS tutorial demo program: Diffusion of a Gaussian hill.
u'= Laplace(u) + f in a square domain
u = u_D on the boundary
u = u_0 at t = 0
u_D = f = 0
The initial condition u_0 is chosen as a Gaussian hill.
"""
from __future__ import print_function
from fenics import *
import time
T = 2.0 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
D0 = 0.1
lambda_1 = -1
# Create mesh and define function space
nx = ny = 30
mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, Constant(0), boundary)
# Define initial value
u_0 = Expression('exp(-a*pow(x[0], 2) - a*pow(x[1], 2))',
degree=2, a=5)
u_n = interpolate(u_0, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
D = Constant(D0)
F = u*v*dx + dt*D*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)
# Create VTK file for saving solution
vtkfile = File('heat_gaussian/solution.pvd')
# Time-stepping
u = Function(V)
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(a == L, u, bc)
# Save to file and plot solution
vtkfile << (u, t)
# Update previous solution
u_n.assign(u)
# Update D_0
D = D0*exp(-lambda_1*t)
D = Constant(D)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
# Update F
F = u*v*dx + dt*D*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
u = Function(V)
```

The initial plot can be displayed perfectly in paraview, but the following images were all blank, not as it should be. Any idea WHY? This is the minimal example of a changing-parameter-simulation using changing diffusion coefficient of heat transpotation.

Community: FEniCS Project

### 1 Answer

0

Solved it with the code below. May it be helpful to others.

```
for n in range(num_steps):
# Update current time
t += dt
# Compute solution
solve(a == L, u, bc)
# Save to file and plot solution
vtkfile << (u, t)
######################################################
# ADDED CODES to ft04 ##
######################################################
# Update D_0
D = D0*np.exp(-lambda_1*t)
D = Constant(D)
u1 = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
# Update F
F = u1*v*dx + dt*D*dot(grad(u1), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)
######################################################
u_n.assign(u)
```

Please login to add an answer/comment or follow this question.