Error: Unable to read vector from file.
*** Reason: Size mis-match between vector in file and input vector.
from __future__ import print_function
from fenics import *
E = 0.14 # Evaporation
T = 5.0 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
M1 = 0.53 # Change in moisture stratification
D = 75000 # Moisture diffusion coefficient
alpha_c = 0.25 # Moisture adjustment time scale
q_s = 75 # Column saturation moisture content
a=0.0
q_c = 60 # Moisture threshold for strong precipitation
alpha = 0.5 #Moisture adjustment time scale
M_s0 = 31400 # Reference value of gross dry stability
Gamma = 4860 # Change in M_s per unit change in u_2
L=1
P_c=0
Q_c=0
F_c=0
F_s=0
M2 = -0.98 # Change in condenstae stratification
# Create mesh and define function space
nx = ny = 256
mesh = RectangleMesh( Point(0, 0), Point(5120, 6400), nx, ny)
# Define function space for velocity
W = VectorFunctionSpace(mesh, 'P', 2)
#Define function space for system of equations
P1 = FiniteElement('P', triangle, 1)
element = MixedElement( [P1, P1] )
V = FunctionSpace ( mesh, element)
# Define test functions
v_1, v_2= TestFunctions(V)
# Define functions for velocity and PDE
w = Function(W)
u = Function(V)
u_n = Function(V)
# Split system functions to access components
u_1, u_2 = split(u)
u_n1, u_n2 = split(u_n)
#Define source terms
f_1 = Constant (0.14)
f_2 = Constant (0.0)
# Define expressions used in variational forms
k = Constant(dt)
M1 = Constant(M1)
D = Constant(D)
alpha_c = Constant(alpha_c)
q_s = Constant(q_s)
q_c = Constant(q_c)
M_s0 = Constant(M_s0)
Gamma = Constant(Gamma)
M2 = Constant(M2)
# Define re-evaporation function
def Re(u_1, u_2):
return ((alpha_c/2)*u_2* (1-(u_1/q_s)))
#Define Latent heat source
def P(u_1):
return ( alpha*(u_1-q_c))
#Define vapor to condensate conversion function
def Q(u_1):
return(2*alpha* (u_1-q_s))
#Define function
def K(u_1, u_2):
return( (L*(P(u_1)-P_c)+L*(Q(u_1)-Q_c)+(F_s-F_c))/M_s(u_2))
def M_s(u_2):
return( M_s0 + Gamma * u_2)
#Define function of large scale precipitation
def P_(u_1, u_2):
return ( (alpha_c/2)*u_2* (u_1/q_s))
# Define Heaviside function H
def H(x):
if ( x>0):
return(1)
else :
return (0)
# Define variational problem
F = ((u_1 -u_n1) / k)*v_1*dx - M1* dot(w, grad(u_1))*v_1*dx \
+ D* dot(grad(u_1), grad(v_1))*dx - Re(u_1, u_2)*v_1*dx \
+ P(u_1)*v_1*dx + Q(u_1)*v_1*dx - M1*K(u_1, u_2)*u_1*v_1*dx \
+ ((u_2 -u_n2) / k)*v_2*dx - M2* dot(w, grad(u_2))*v_2*dx \
+ D* dot(grad(u_2), grad(v_2))*dx - Q(u_1)*v_1*dx \
+ P_(u_1, u_2)*v_2*dx - M2*K(u_1, u_2)*u_2*v_2*dx \
- f_1*v_1*dx - f_2*v_2*dx
# Create time series for reading velocity data
timeseries_w = TimeSeries( 'navier_stokes_cylinder/velocity_series')
#timeseries_w = TimeSeries('navier_stokes_cylinder/velocity_series')
#Create VTK files for visualization output
vtkfile_u_1 = File( 'reaction_system/u_1.pvd')
vtkfile_u_2 = File( 'reaction_system/u_2.pvd')
# Create progress bar
progress = Progress('Time-stepping')
# Time-stepping
t = 0
for n in range(num_steps):
t += dt
timeseries_w.retrieve(w.vector(), t)
solve(F == 0, u)
_u_1, _u_2 = u.split()
vtkfile_u_1 << (_u_1, t)
vtkfile_u_2 << (_u_2, t)
u_n.assign(u)
steps = 500
p = Progress("Looping", steps)
for i in range(steps):
set_log_level(LogLevel.PROGRESS)
p+=1
set_log_level(LogLevel.ERROR)
#progress.update(t / T)
#interactive()
# Plot solution
import matplotlib.pyplot as plt
plot(mesh)
plot(_u_1)
plot(_u_2)
plt.show()