3 component advection diffusion reaction solver added to darcy advect…

…ion class and a new script LVL_3comp added

src/modules/mupopp.py

@@ -848,7 +848,117 @@ def advection_diffusion_two_component(self,W,mesh,sol_prev,dt,f1,K=1.0,zh=Consta
         term_SUPG = tau_SUPG*dot(u, grad(vc))*r*dx
         F += term_SUPG       
         return lhs(F), rhs(F)
-
+    def advection_diffusion_three_component(self,W,mesh,sol_prev,dt,f1,K=1.0,zh=Constant((0.0,1.0)),SUPG=1,gam_rho=-1.0,rho_0=1.0):
+        """     
+        This function returns the bilinear form for the system of PDES governing
+        an advection-reaction-two-component flow, the governing
+        PDEs for Darcy flow are:        
+        div(u) = 0                                                          (1)
+        phi*u = -k*(grad(p)-rho*zhat)                                      (2)
+        dc0/dt + dot(u,grad(c0)) = div(grad(c0))/Pe - Da*c0*c1**2/phi + beta*f (3)
+        and  
+        dc1/dt = - 2*Da*c0*c1**2/phi                                             (4)
+        where 
+        drho = 1 + gam_rho c0                                                       (5)
+        where c0 and c1 are concentrations of the reactants in the liquid
+        and solid, u is the fluid velocity, p is pressure, k is permeability
+        drho=difference between liquid and solid densities, zhat is a unit
+        vecotr in vertically upward direction, Pe is Peclet number, Da is
+        the Dahmkoler number, beta is source strength, f is a function
+        for lateral variations in source of c0, and phi is the constant porosity
+        Input:
+             W        : Mixed function space containing velocity, pressure and two
+                        concentrations
+             mesh     : Fenics mesh on which W is defined
+             sol_prev : All unknowns from the previous time step
+             dt       : time step
+             f1       : Spatially variable scalar function for the source term
+             K        : Constant permeability
+             zh       : Unit vector in the vertical direction
+             SUPG     : A counter for chosing between different SUPG terms
+             gam_rho  : The coefficient of concentration in the expression for drho
+             rho_0    : The intercept in the expression for drho
+        Output:
+             lhs(F)   : Left hand side of the combined bilinear form
+             rhs(F)   : Right hand side of the bilinear formulation
+        The bilinear form uses Crank-Nicholson time stepping and gives 
+        several options for SUPG stabilization. We recommend using Sendur 2018 formulation
+        for the SUPG term.
+        The bilinear form, combined together, becomes
+        phi*dot(u,v)-k*p*div(v)+div(u)*q + eta*(c0n-c0nm1)+eta*dot(u,grad(c_mid))*dt
+                 +dt*dot(grad(c_mid),grad(eta))/Pe+(c1_n-c1_nm1)*omega
+                 = K*drho*dot(zhat,v)-Da*c0*c1*dt*eta/phi 
+                   - beta*f1*dt*eta + Da*c0*c1*dt*omega/(1-phi)
+        In this code c0 is uc, c1 is cc, eta is vc, omega is qc
+        We recommend using this function preferentially over the other functions as
+        this combines all four equations into one bilinear form. If the density contrast
+        is constant, just replace drho with a constant value.
+        """
+	h = CellDiameter(mesh)
+	# TrialFunctions and TestFunctions
+        U = TrialFunction(W)
+        (v, q, vc, qc,qc1) = TestFunctions(W)
+        u, p, uc, cc,cc1   = split(U)
+        zhat=zh
+        # Density contrast as a function of concentration
+        deltarho=rho_0-gam_rho*uc
+
+        # Define the variational form
+        F = (inner(self.phi*u,v) - K*div(v)*p+div(u)*q)*dx\
+            - K*deltarho*inner(v,zhat)*dx
+        # uc and cc are the trial functions for the next time step
+        # uc for comp cc and d comp1 
+        # u0 (component 0), c0(component 1), and c1(component 3)
+        # are known values from the previous time step
+        u,p,u0 ,c0,c1 = split(sol_prev)
+        # Mid-point solution for comp 0
+        u_mid = 0.5*(u0 + uc)
+        # First order reaction term
+        # For chemical reaction
+        # MgCO3+2Fe = C + 3 Fe(2/3)Mg(1/3)O
+        # Assuming forward reaction controls the rate
+        # Total molar mass is calculated from
+        # MgCO3 = 24.0+12.0+3.0*16.0
+        # 2Fe = 2.0*56
+        # C   = 12.0
+        # 3 Fe(2/3)Mg(1/3)O = 3.0*(0.67*56.0+0.33*24.0+16.0)
+        # M = MgCO3+2*Fe+C+3*Oxide = 392.32
+        Total_molar_mass = 392.32
+        Fe_frac = 2.0*56.0/Total_molar_mass
+        carbonate_frac = 84.0/Total_molar_mass
+        C_frac = 12.0/Total_molar_mass
+        f21 = Fe_frac*self.Da*u0*c0*c0 #reaction rate for Fe
+        f11 = carbonate_frac*self.Da*u0*c0*c0 #reaction rate for carbonate
+        f31 = C_frac*self.Da*u0*c0*c0
+
+	F += vc*(uc - u0)*dx + dt*(vc*dot(u, grad(u_mid))*dx\
+                + dot(grad(vc), grad(u_mid)/self.Pe)*dx) \
+		+ dt*f11/self.phi*vc*dx  - self.beta*dt*f1*vc*dx \
+                + qc*(cc - c0)*dx + dt*f21/(1-self.phi)*qc*dx \
+                + qc1*(cc1 - c1)*dx - dt*(f31/(1-self.phi))*qc1*dx
+        # Residual
+        h = CellDiameter(mesh)
+        r = uc - u0 + dt*(dot(u, grad(u_mid)) - div(grad(u_mid))/self.Pe+f11/self.phi)\
+            - self.beta*dt*f1 + cc-c0 + dt*f21/(1.0-self.phi) + cc1-c1 - dt*f31/(1.0-self.phi)
+        # Add SUPG stabilisation terms
+        # Default is Sendur 2018, a modification of Codina, 1997 also works
+        vnorm = sqrt(dot(u, u))
+
+        if SUPG==1:
+            #Sendur 2018
+            tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h)
+        elif SUPG==2:
+            #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*np.max(c0)/self.phi)
+        else:
+            alpha_SUPG=self.Pe*vnorm*h/2.0
+            #Brookes and Hughes
+            coth = (np.e**(2.0*alpha_SUPG)+1.0)/(np.e**(2.0*alpha_SUPG)-1.0)
+            tau_SUPG=0.5*h*(coth-1.0/alpha_SUPG)/vnorm
+        term_SUPG = tau_SUPG*dot(u, grad(vc))*r*dx
+        F += term_SUPG       
+        return lhs(F), rhs(F)
 
 ###################################################################
 ### Advection diffusion equation in Stokes flow

src/run/LVL_3comp/LVL_3comp.py

@@ -0,0 +1,246 @@
+###################################################
+## This is a part of Mupopp
+## Copyright Saswata Hier-Majumder, July, 2016
+## Modified on August 2018
+## This program solves an advection diffusion problem
+## with Darcy flow, using Dirichlet boundary conditions
+## for velocity and concentration
+## Modified by Joe Sun, February 2019
+####################################################
+
+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 *
+
+
+#####################################################
+parameters["std_out_all_processes"]=False
+#set_log_level(30) #Critical info only
+####################################
+
+# Parameters for initializing the object
+Da0  = 100.0
+Pe0  = 100.0
+alpha0= 0.01   # beta = alpha0/phi
+phi0 = 0.01
+beta=alpha0/phi0
+Fe=0.01
+cfl0 = 0.01
+# Parameters for iteration
+T0 = 2.0
+dt0 = 0.1
+out_freq0 = 1
+
+# Parameters for mesh
+mesh_density = 40
+
+# Output files for quick visualisation
+
+file_name       =  "Da_%3.2f_Pe_%.1E_beta_%3.2f_Fe_%3.2f"%(Da0,Pe0,beta,Fe)
+output_dir     =  "output/"
+
+extension      = "pvd"   # "xdmf" or "pvd"
+
+velocity_out   = File(output_dir + file_name + "_velocity." + extension, "compressed")
+pressure_out   = File(output_dir + file_name + "_pressure." + extension, "compressed")
+c0_out         = File(output_dir + file_name + "_concentration0." + extension, "compressed")
+c1_out         = File(output_dir + file_name + "_concentration1." + extension, "compressed")
+c2_out         = File(output_dir + file_name + "_concentration2." + extension, "compressed")
+initial_c0_out = File(output_dir + file_name + "_initial_c0." + extension, "compressed")
+initial_c1_out = File(output_dir + file_name + "_initial_c1." + extension, "compressed")
+
+# Output parameters
+def output_write(mesh_density,Da,phi,Pe,alpha,cfl,fname= output_dir + "/a_parameters.out"):
+    """This function saves the output of iterations"""
+    file=open(fname,"a")
+    file.write("####################################")
+    file.write("\n")
+    file.write("Mesh density:  %g\n" %mesh_density)
+    file.write("Da:  %g\n" %Da0)
+    file.write("phi:  %g\n" %phi0)
+    file.write("Pe:  %g\n" %Pe0)
+    file.write("alpha:  %g\n" %alpha0)
+    file.write("cfl:  %g\n" %cfl0)
+    file.write("####################################")
+    file.close
+
+output_write(mesh_density,Da0,phi0,Pe0,alpha0,cfl0)
+
+#Define function for source term in Governing equations
+class SourceTerm(Expression):
+    """ Creates an expression for the source term
+    in the advection reaction equations.
+    The source term consists of a series of
+    sine waves.
+    """
+    def __init__(self, mesh,element):
+        self.mesh = mesh
+        self.element=element
+    def eval(self, values, x):
+        g1=x[0]*0.0
+        for ii in range(0,20):
+            g1+=0.1*np.abs(np.sin(ii*x[0]*np.pi))
+
+        g = (1.0-tanh(x[1]/0.01))*g1
+        values[0] = g
+    def value_shape(self):
+        return (1,)
+
+# Define function for BC
+class BoundarySource(Expression):
+    def __init__(self, mesh,element):
+        self.mesh = mesh
+        self.element=element
+    def eval_cell(self, values, x, ufl_cell):
+        cell = Cell(self.mesh, ufl_cell.index)
+        n = cell.normal(ufl_cell.local_facet)
+        g1=x[0]*0.0
+        for ii in range(0,20):
+            g1+=0.01*np.abs(np.sin(ii*x[0]*np.pi))
+        g = -g1
+        values[0] = g*n[0]
+        values[1] = g*n[1]
+    def value_shape(self):
+        return (2,)
+
+
+############################
+## 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))
+mesh   = generate_mesh(domain,mesh_density)
+#mesh = RectangleMesh(Point(xmin, ymin), Point(xmax, ymax), 100, 50)
+#####################
+
+# Mark facets
+## Define subdomain for flux calculation
+class Plane(SubDomain):
+  def inside(self, x, on_boundary):
+    return x[1] > ymax - DOLFIN_EPS
+facets = FacetFunction("size_t", mesh)
+Plane().mark(facets, 1)
+ds = Measure("ds")[facets]
+n = FacetNormal(mesh)
+####################
+# 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()            
+
+############################
+## Darcy velocity
+############################
+# Create FunctionSpaces
+# Velocity
+V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
+# Pressure
+Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
+#Concentration
+Qc = FiniteElement("Lagrange",mesh.ufl_cell(), 1)
+# Make a mixed space
+W = dolfin.FunctionSpace(mesh, MixedElement([V,Q,Qc,Qc,Qc]), constrained_domain=pbc)
+X  = FunctionSpace(mesh,"CG",1, constrained_domain=pbc)
+
+
+# Define boundary conditions
+G=BoundarySource(mesh,element=V)
+#bc1 = DirichletBC(W.sub(0), G, bottom)
+bc1 = DirichletBC(W.sub(0), Constant((0.0,0.1)), bottom)
+bc  = [bc1]
+
+###########################
+## Create an object
+###########################
+darcy = DarcyAdvection(Da=Da0,phi=phi0,Pe=Pe0,alpha=alpha0,cfl=cfl0)
+
+# Define initial conditions
+sol_0 = Function(W)
+temp2 = Function(X)
+c01_temp=Expression("0.01",element=Q) #Fe
+temp2.interpolate(c01_temp)
+c01 = temp2
+assign(sol_0.sub(3), c01)
+
+#c0_initial,c1_initial=c0.split()
+#initial_c0_out << c0_initial
+#initial_c1_out << c1_initial
+
+###########################
+## Solve for Darcy velocity
+###########################
+sol = Function(W)
+# Parameters for iteration
+T = T0
+dt = dt0
+t = dt
+flux=np.array([])
+carbon=np.array([])
+time_array=np.array([])
+i = 1
+out_freq = out_freq0
+S=SourceTerm(mesh,element=Qc)
+
+while t - T < DOLFIN_EPS:
+    # Update the concentration of component 0
+    a,L = darcy.advection_diffusion_three_component(W,mesh,sol_0,dt,f1=S,K=0.1)
+    solve(a==L,sol,bc)
+    sol_0 = sol
+    u0,p0,c00,c01,c02 = sol.split()
+    if i % out_freq == 0:
+	u0.rename("velocity","")
+	velocity_out << u0
+	p0.rename("pressure","")
+	pressure_out << p0
+        c00.rename("[CO3]","")
+        c0_out << c00
+        c01.rename("[Fe]","")
+        c1_out << c01
+        c02.rename("[C]","")
+        c2_out << c02
+        ## Calculate flux
+        flux1 = assemble(c00*phi0*dot(u0, n)*ds(1))
+        flux= np.append(flux,flux1)
+        time_array=np.append(time_array,t)
+        #calculate mass of carbon deposited
+        mass1=assemble(c02*(1.0-phi0)*dx)
+        carbon = np.append(carbon,mass1)
+        #print "flux 1: ", flux_1
+    # Move to next interval and adjust boundary condition
+    info("time t =%g\n" %t)
+    info("iteration =%g\n" %i)
+    #print 'iteration',i
+    t += dt
+    i += 1
+flux_file=output_dir + file_name + "_flux.csv"
+np.savetxt(flux_file,(time_array,flux),delimiter=',')
+mass_file=output_dir + file_name + "_mass.csv"
+np.savetxt(mass_file,(time_array,carbon),delimiter=',')