Medium modification in SUPG stabilization term from Cordina

src/modules/mupopp.py

@@ -656,12 +656,12 @@ def advection_diffusion_two_component_nonadaptive(self,Q, c_prev, velocity, dt,m
         vnorm = sqrt(dot(velocity, velocity))
         #F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx
         alpha_SUPG=self.Pe*vnorm*h/2.0
-        term1_SUPG=0.5*h*(np.arctanh(alpha_SUPG)-1.0/alpha_SUPG)/vnorm
+        term1_SUPG=0.5*h*(1.0/np.tanh(alpha_SUPG)-1.0/alpha_SUPG)/vnorm
         term_SUPG = term1_SUPG*dot(velocity, grad(v))*r*dx
         a = 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
-        L =     dt*f*v*dx + dt*f*q*self.phi*dx
+        L =     -dt*f*v*dx - dt*f*q*self.phi*dx
 
         #return lhs(F), rhs(F)
         return a, L
@@ -713,16 +713,38 @@ def advection_diffusion_two_component(self,Q, c_prev, velocity,mesh):
         #####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
         # Residual
         r = u - u0 + dt*(dot(velocity, grad(u_mid)) - div(grad(u_mid))/self.Pe+f)+c-c0+dt*f
-        # Galerkin variational problem
+        # Add SUPG stabilisation terms
+        vnorm = sqrt(dot(velocity, velocity))
+        #F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx
+
+        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)
+        term1_SUPG=0.5*h*(coth-1.0/alpha_SUPG)/vnorm
+        #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)
+        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)+dt*f*v*dx\
-                + q*(c - c0)*dx + dt*f*q*self.phi*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 
+
+        # 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
 
         # Add SUPG stabilisation terms
-        vnorm = sqrt(dot(velocity, velocity))
-        F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx 
+        #vnorm = sqrt(dot(velocity, velocity))
+        #F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx 
 
         return lhs(F), rhs(F)
 

src/run/LVL_RII/LVL_RII.py

@@ -106,7 +106,7 @@ def map(self, x, y):
 bc=[bc1,bc2]
 
 ##############Create an object
-darcy=DarcyAdvection(Da=10.0,phi=0.01,Pe=1.0e4,cfl=0.1)
+darcy=DarcyAdvection(Da=10.0,phi=0.01,Pe=1.0e6,cfl=0.1)
 
 ########################
 ## Solve for Darcy velocity
@@ -143,13 +143,13 @@ def map(self, x, y):
 c00=temp
 c00.rename("[CO3]","")
 initial_c0_out << c00
-temp.interpolate(Expression("0.9",degree=1))
+temp.interpolate(Expression("0.1",degree=1))
 c01=temp
 # Write the initial concentrations into files
 
 
 # Parameters
-T = 5.0
+T = 50.0
 darcy.dt = 0.10
 t = darcy.dt
 
@@ -160,7 +160,7 @@ def map(self, x, y):
 #bc2 = DirichletBC(Q1.sub(0), Constant(0.5),bottom)
 # Set up boundary conditions for component 1
 #bc3 = DirichletBC(Q1.sub(1), Constant(0.9),top)
-bc4 = DirichletBC(Q1.sub(1), Constant(0.9),bottom)
+bc4 = DirichletBC(Q1.sub(1), Constant(0.1),bottom)
 
 #bc_c=[bc1,bc2,bc3,bc4]
 bc_c=[bc2,bc4]
@@ -199,4 +199,4 @@ def map(self, x, y):
 plt.loglog(time_array,dt_array,'or')
 plt.xlabel('time')
 plt.ylabel('dt')
-plt.show()
+#plt.show()