[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)

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
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
    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.

Similar posts:
Search »