From e33c6f60f2f3410f340b366b6c96f9ca35fdda29 Mon Sep 17 00:00:00 2001 From: sashgeophysics Date: Tue, 2 Oct 2018 15:46:39 +0100 Subject: [PATCH] Medium modification in SUPG stabilization term from Cordina --- src/modules/mupopp.py | 36 +++++++++++++++++++++++++++++------- src/run/LVL_RII/LVL_RII.py | 10 +++++----- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/src/modules/mupopp.py b/src/modules/mupopp.py index eb96255..bfbf5d1 100644 --- a/src/modules/mupopp.py +++ b/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) diff --git a/src/run/LVL_RII/LVL_RII.py b/src/run/LVL_RII/LVL_RII.py index 8d8e5cd..b8b4c58 100644 --- a/src/run/LVL_RII/LVL_RII.py +++ b/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()