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

157
views
1
3 months ago by
# 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)

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

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

0
3 months ago by

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