change in solution curve with change in time step size
7 weeks ago by
I posted a question yesterday that was viewed a fair amount, but was never answered. Just in case it was a poorly phrased question, I'll rewrite it.
I've adjusted the advection-diffusion-reaction example code (included below) to monitor the flow of molecules out of the system. I would expect the observed curve to converge to a single curve with as the size of the time step is decreased, but it seems to continually change (it doesn't stop shrinking). I'm not sure if I'm not understanding something fundamental or if I've written something incorrectly in the code. Thanks again!
from fenics import * import numpy as np import matplotlib.pyplot as plt import random import sys import time import warnings fig2, ax2 = plt.subplots() ax2.set_xlabel('$t$') ax2.set_ylabel('$C_i@exit$') time_step = [100,500,1000,5000,10000,50000,100000,500000,1000000] for num in time_step: T = 5 num_steps = num dt = T / num_steps #k = Constant(dt) set_log_active(False) warnings.filterwarnings("ignore", category=DeprecationWarning) #T = 5.0 # final time #num_steps = 500 # number of time steps #dt = T / num_steps # time step size eps = 0.01 # diffusion coefficient K = 10.0 # reaction rate mesh_size = 100 # Read mesh from file #mesh = Mesh('navier_stokes_cylinder/cylinder.xml.gz') mesh = UnitIntervalMesh(mesh_size) #W = VectorFunctionSpace(mesh,'P',2) #element = MixedElement([P1]) V = FunctionSpace(mesh,'P',1) ## Define function space for velocity #W = VectorFunctionSpace(mesh, 'P', 2) # Define function space for system of concentrations #P1 = FiniteElement('P', triangle, 1) #element = MixedElement([P1]) #V = FunctionSpace(mesh, element) # Define test functions v_1 = TestFunction(V) # Define functions for velocity and concentrations #w = Function(W) u = Function(V) u_n = Function(V) # Split system functions to access components #u_1 = split(u) #u_n1 = split(u_n) # Define expressions used in variational forms k = Constant(dt) eps = Constant(eps) #dx = Measure("dx")[domains] # Define variational problem F = ((u - u_n) / k)*v_1*dx+ dot(grad(u), grad(v_1))*dx # Create time series for reading velocity data #F = ((u_1 - u_n1) / k)*v_1*dx + dot(w, grad(u_1))*v_1*dx #F = ((u_1 - u_n1) / k)*v_1*dx(0) + eps*dot(grad(u_1), grad(v_1))*dx(0) # Create progress bar progress = Progress('Time-stepping') set_log_level(PROGRESS) # Time-stepping t = 0 tol = 1e-14 def boundary_L(x, on_boundary): return on_boundary and near(x,0,tol) def boundary_R(x, on_boundary): return on_boundary and near(x,1,tol) bc_R_1 = DirichletBC(V,Constant(0),boundary_R) #bc_R_2 = DirichletBC(V,Constant(0),boundary_R) conVtime =  conVtime_2 =  timing =  timing.append(t) conVtime.append(0) for n in range(num_steps): # Update current time t += dt print("TIME") print(t) # Read velocity from file if t > dt: bcs = [bc_R_1] else: u_L = (1./(2*dt))*(1+np.sin((np.pi/dt)*(t-dt/2))) #u_L = 10000 bc_L_1 = DirichletBC(V,Constant((u_L)),boundary_L) bcs = [bc_L_1, bc_R_1] print(u_L) # Solve variational problem for time step solve(F == 0, u, bcs) # Save solution to file (VTK) timing.append(t) conVtime.append(u.vector().get_local()) # Update previous solution u_n.assign(u) # Update progress bar progress.update(t / T) #sys.exit() ax2.plot(timing,conVtime, ls='--',color='r', alpha=0.7) plt.show()
Community: FEniCS Project
7 weeks ago by
Probably your boundary conditions are not such that you should have expected a convergence (if there is any). Have better look at u_L, it depends on dt as well! Moreover, I am not sure if the pde that you tend to solve has to converge in the way that you want to see. I would try for the steady state case (if the velocity is slow then you will find the solution).
Please login to add an answer/comment or follow this question.