change in solution curve with change in time step size


69
views
0
7 weeks ago by
Hi again,

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],0,tol)

  def boundary_R(x, on_boundary):
        return on_boundary and near(x[0],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()[1])

      # 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

1 Answer


1
7 weeks ago by
Emek  
It is really hard to understand your code with lots of comment outs and different ideas. Please consider that the reader does not know the details at all.

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

Best, Emek
Hi Emek,

Thank you for the response. Yes, I really should have taken the time to clean up the code and will do so in the future. Sorry about that!

I'll try these suggestions out and see if I make any progress.

Best,
Adam
written 7 weeks ago by Adam Yonge  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »