advection test added and mupopp modified

src/modules/mupopp.py

@@ -564,14 +564,15 @@ class DarcyAdvection():
     
     """
 
-    def __init__(self,Pe=100,Da=10.0,phi=0.01,cfl=1.0e-2,dt=1.0e-2):
+    def __init__(self,Pe=100,Da=10.0,phi=0.01,alpha=0.005,cfl=1.0e-2,dt=1.0e-2):
         """Initiates the class. Inherits nondimensional numbers
         from compaction, only need to add Pe"""
         self.Pe=Pe
         self.Da=Da
         self.phi=phi
         self.cfl=cfl
         self.dt=dt
+        self.alpha=alpha
     def darcy_bilinear(self,W,mesh,K=0.1,zh=Constant((0.0,1.0))):
         """            
         """
@@ -672,7 +673,8 @@ def advection_diffusion_two_component(self,Q, c_prev, velocity,mesh):
         calculated outside this function after solving the bilinear form from
         this step"""
 
-        f  = Expression("0.0",degree=1)        
+        f1  = Expression("0.005*(1.0-tanh(x[1]/0.2))*sin(2.0*x[0]*3.14)",degree=1)
+        f2  = Expression("0.001",degree=1)
         h = CellDiameter(mesh)
         # Parameters
         U = TrialFunction(Q)
@@ -691,6 +693,13 @@ def advection_diffusion_two_component(self,Q, c_prev, velocity,mesh):
 
         # First order reaction term
         f = self.Da*u0*c0
+        # the source term for c0 f1=[0,1]
+	f1 = Expression('(1.0-tanh(x[1]/0.01))*( (1.0+sin(1.0*x[0]*3.14)) \
+		+ (1.0+sin(2.0*x[0]*3.14)) + (1.0+sin(3.0*x[0]*3.14)) \
+		+ (1.0+sin(4.0*x[0]*3.14)) + (1.0+sin(5.0*x[0]*3.14)) \
+		+ (1.0+sin(6.0*x[0]*3.14)) + (1.0+sin(7.0*x[0]*3.14)) \
+		+ (1.0+sin(8.0*x[0]*3.14)) + (1.0+sin(9.0*x[0]*3.14)) \
+		+ (1.0+sin(10.0*x[0]*3.14)) )/20.0',degree=1)
         ###############################################
         # Adaptive time-stepping added September 2018
         # Right hand side of advection equation
@@ -705,20 +714,21 @@ def advection_diffusion_two_component(self,Q, c_prev, velocity,mesh):
         ad_max=advect_rhs.vector().max()
         #Compute dt for this time step such that
         #dt*ad_max<=CFL
-        if np.abs(ad_max)>self.cfl:
-            self.dt=self.cfl/ad_max
-
+        #if np.abs(ad_max)>self.cfl:
+        #    self.dt=self.cfl/ad_max      
         dt=self.dt
-        #print 'dt',dt,'ad_max',ad_max
         #####End adaptive time stepping
         ##################################################
 
         # Galerkin variational problem
         #F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx\
-        #        + dot(grad(v), grad(u_mid)/self.Pe)*dx)+dt*f*v*dx\
-        #        + q*(c - c0)*dx + dt*f*q*self.phi*dx
+        #        + dot(grad(v), grad(u_mid)/self.Pe)*dx)+(dt*f*v/self.phi)*dx\
+        #        + q*(c - c0)*dx + dt*f*q*dx
         # Residual
-        r = u - u0 + dt*(dot(velocity, grad(u_mid)) - div(grad(u_mid))/self.Pe+f)+c-c0+dt*f
+        #r = u - u0 + dt*(dot(velocity, grad(u_mid)) - div(grad(u_mid))/self.Pe+self.Da*u0*c0)\
+        #    +c-c0+dt*self.Da*u0*c0/self.phi
+        r = u - u0 + dt*(dot(velocity, grad(u_mid)) - div(grad(u_mid))/self.Pe+f/self.phi)\
+            +c-c0+dt*f/(1.0-self.phi)
         # Add SUPG stabilisation terms
         vnorm = sqrt(dot(velocity, velocity))
         #F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx
@@ -730,12 +740,13 @@ def advection_diffusion_two_component(self,Q, c_prev, velocity,mesh):
         #Sendur 2018
         tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h)
         #Codina 1997 eq. 114
-        tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h+self.Da)
+        #tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h+self.Da)
+        tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h+self.Da*np.max(c0)/self.phi)
         term_SUPG = tau_SUPG*dot(velocity, grad(v))*r*dx
-        F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx\
-                + dot(grad(v), grad(u_mid)/self.Pe)*dx)\
-                + q*(c - c0)*dx  + term_SUPG
-        F +=     -dt*f*v*dx - dt*f*q*self.phi*dx 
+        #F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx\
+        #        + dot(grad(v), grad(u_mid)/self.Pe)*dx)\
+        #        + q*(c - c0)*dx  + term_SUPG
+        #F +=    f1*v*dx+f2*q*dx -(dt*self.Da*u0*c0*v/self.phi)*dx - dt*self.Da*u0*c0*q*dx 
 
         # Galerkin variational problem
         #F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx\
@@ -745,7 +756,12 @@ def advection_diffusion_two_component(self,Q, c_prev, velocity,mesh):
         # Add SUPG stabilisation terms
         #vnorm = sqrt(dot(velocity, velocity))
         #F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx 
-
+        F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx \
+                                + dot(grad(v), grad(u_mid)/self.Pe)*dx) \
+                                + dt*f/self.phi*v*dx  - self.alpha/self.phi*dt*f1*v*dx\
+                                + q*(c - c0)*dx + dt*f/(1-self.phi)*q*dx
+        #5.f/self.phi 6.dt*f/(1-self.phi) 9.f1 
+        F += term_SUPG
         return lhs(F), rhs(F)
 
 

src/run/LVL_RII/advection_test.py

@@ -0,0 +1,189 @@
+###################################################
+## This is a part of Mupopp
+## Copyright Saswata Hier-Majumder, July, 2016
+## Modified on December 2018
+## This program solves an advection diffusion problem
+## with Darcy flow, using Dirichlet boundary conditions
+## for velocity and concentration
+####################################################
+
+from fenics import *
+from mshr import*
+import numpy, scipy, sys, math
+#import matplotlib.pyplot as plt
+#Add the path to Mupopp module to the code
+sys.path.insert(0, '../../modules/')
+#from mupopp import *
+
+# muppop_Joe contains the newest equations
+from mupopp import *
+
+class BoundaryConcentration(Expression):
+    def __init__(self, mesh, element):
+        self.mesh = mesh
+        self.element=element
+    def eval_cell(self, values, x, ufl_cell):
+        f = 0.05*(1.0-tanh(0.0/0.2))*(1.0+sin(2.0*x[0]*3.14))/2.0
+        values[0] = f       
+    def value_shape(self):
+        return (1,)
+#####################################################
+parameters["std_out_all_processes"]=False
+#set_log_level(30) #Critical info only
+# Parameters for initializing the object
+Da0  = 1.0
+phi0 = 0.01
+Pe0  = 1.0e4
+alpha0= 0.01   # beta = alpha0/phi
+cfl0 = 0.1
+v_temp=Expression(("0.0","0.1"),degree=2)
+c01_temp=Expression("0.01",degree=1)
+
+# Parameters for iteration
+T0 = 0.1
+dt0 = 1.0e-4
+out_freq0 = 5
+####################################
+# Output files for quick visualisation
+output_dir     = "output/"
+extension      = "pvd"   # "xdmf" or "pvd"
+
+velocity_out   = File(output_dir + "velocity." + extension, "compressed")
+pressure_out   = File(output_dir + "pressure." + extension, "compressed")
+c0_out          = File(output_dir + "concentration0." + extension, "compressed")
+c1_out          = File(output_dir + "concentration1." + extension, "compressed")
+initial_c0_out   = File(output_dir + "initial_c0." + extension, "compressed")
+initial_c1_out   = File(output_dir + "initial_c1." + extension, "compressed")
+############################
+## Numerical solution
+############################
+
+# Define the mesh
+xmin = 0.0
+xmax = 4.0
+ymin = 0.0
+ymax = 1.0
+domain = Rectangle(Point(xmin,ymin),Point(xmax,ymax))
+# Parameters for mesh
+mesh_density = 60
+mesh   = generate_mesh(domain,mesh_density)
+#mesh = RectangleMesh(Point(xmin, ymin), Point(xmax, ymax), 100, 50)
+
+# Define essential boundary
+def top_bottom(x):
+    return x[1] < DOLFIN_EPS or x[1] > ymax - DOLFIN_EPS
+def bottom(x):
+    return x[1] < DOLFIN_EPS
+def top(x):
+    return x[1] > ymax - DOLFIN_EPS
+
+# Sub domain for Periodic boundary condition
+class PeriodicBoundary(SubDomain):
+    # Bottom boundary is "target domain" G
+    def inside(self, x, on_boundary):
+        return bool(near(x[0], xmin) and on_boundary)   
+    # Map right boundary (H) to left boundary (G)
+    def map(self, x, y):
+        if near(x[0], xmax):
+            y[0] = x[0] - xmax
+            y[1] = x[1]            
+        else:
+            y[0] = -1000
+            y[1] = -1000
+# Create periodic boundary condition
+pbc = PeriodicBoundary()            
+####################################################
+# Create FunctionSpaces
+Qc = FiniteElement("Lagrange",mesh.ufl_cell(), 1)
+QM = dolfin.FunctionSpace(mesh, MixedElement([Qc,Qc]))#,constrained_domain=pbc)
+X  = FunctionSpace(mesh,"CG",1)#, constrained_domain=pbc)
+Vc = VectorFunctionSpace(mesh,"CG",2)#, constrained_domain=pbc)
+
+# Define the constant velocity
+vel = Function(Vc)
+vtemp = Function(Vc)
+vtemp.interpolate(v_temp)
+vel = vtemp
+
+###########################
+## Create an object
+###########################
+darcy = DarcyAdvection(Da=Da0,phi=phi0,Pe=Pe0,alpha=alpha0,cfl=cfl0)
+
+############################
+## Define initial conditions 
+############################
+# Create functions for initial conditions
+c0 = Function(QM)
+
+# Functions for components
+temp1 = Function(X)
+temp1.interpolate(Expression("0.5+0.5*(1.0-tanh(x[1]/0.2))*sin(2.0*x[0]*3.14)",degree=2))
+c00 = temp1
+temp2 = Function(X)
+temp2.interpolate(c01_temp)
+c01 = temp2
+
+# Assign the components
+#assign(c0.sub(0), c00)
+assign(c0.sub(1), c01)
+
+#assign(c0, [c01,c01])
+
+#c_initial = InitialConcentration(element=MixedElement([Qc,Qc]))
+#c0.interpolate(c_initial)
+
+c0_initial,c1_initial=c0.split()
+initial_c0_out << c0_initial
+initial_c1_out << c1_initial
+
+############################
+## Define boundary conditions
+############################
+# Set up boundary conditions for component 0
+c0_bottom = BoundaryConcentration(mesh,element=Qc)
+bc1 = DirichletBC(QM.sub(0), Constant(0.0),top)
+bc2 = DirichletBC(QM.sub(0), c0_bottom,bottom)
+
+# Set up boundary conditions for component 1
+bc3 = DirichletBC(QM.sub(1), Constant(0.1),top)
+bc4 = DirichletBC(QM.sub(1), Constant(0.1),bottom)
+
+#bc_c = [bc1,bc2,bc3,bc4]
+bc_c = [bc2]
+
+###########################
+## Solve for Concentrations
+###########################
+sol_c = Function(QM)
+# Parameters for iteration
+T = T0
+darcy.dt = dt0
+t = darcy.dt
+
+i = 1
+out_freq = out_freq0
+# Save the dt and time values in an array
+dt_array = []
+time_array = []
+while t - T < DOLFIN_EPS:
+
+    # Update the concentration of component 0
+    ac,Lc = darcy.advection_diffusion_two_component(QM,c0,vel,mesh)
+    solve(ac==Lc,sol_c)#,bc_c) # tranfer the DirichletBC into a source term
+    c0 = sol_c
+    #c0.assign(sol_c) # 'assign' is verified to be equal to '=' here
+    c00,c01 = sol_c.split()
+    if i % out_freq == 0:
+        c00.rename("[CO3]","")
+        c0_out << c00
+        c01.rename("[Fe]","")
+        c1_out << c01
+    # Move to next interval and adjust boundary condition
+    dt_array.append(darcy.dt)
+    time_array.append(t)
+    info("time t =%g\n" %t)
+    info("iteration =%g\n" %i)
+    #print 'iteration',i
+    t += darcy.dt
+    i += 1