I tried to find the minimum between the u_1 and the constant q_s. The u_1 defined below in the variational problem.
from __future__ import print_function
from fenics import *
E = 0.14 # Evaporation
T = 5.0 # final time
num_steps = 5 # 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-(min(u_1, q_s)/q_s)))
def Min(u_1, q_s):
return min(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_Ahmed/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)
I defined one function Re(u_1, u_2) and its definition is
def Re(u_1, u_2):
return ((alpha_c/2)u_2 (1-(min(u_1, q_s)/q_s)))
The minimum function does not work for u_1 and q_s. I believe u_1 is in vector form but q_s is a constant value. Is there any way to fix this issue?
The following error occurred.
File “exp.py”, line 101, in
- D* dot(grad(u_1), grad(v_1))*dx - Re(u_1, u_2)v_1dx
File “exp.py”, line 67, in Re
return ((alpha_c/2)u_2 (1-(Min(u_1, q_s)/q_s)))
File “exp.py”, line 70, in Min
return min(u_1,q_s)
File “/Users/kapilchawla/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/conditional.py”, line 46, in bool
error(“UFL conditions cannot be evaluated as bool in a Python context.”)
File “/Users/kapilchawla/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: UFL conditions cannot be evaluated as bool in a Python context.