Add adaptive vscale to stabilization

src/ambit_fe/fluid/fluid_main.py

@@ -274,6 +274,10 @@ def __init__(self, pbase, io_params, time_params, fem_params, constitutive_model
         else:
             self.pre = False
 
+        self.re = fem.Function(self.V_scalar)
+        self.re_ktilde = fem.Function(self.V_v)
+        self.re_c = fem.Function(self.V_v)
+
         # collect references to pressure vectors
         self.pvecs_, self.pvecs_old_ = [], []
         if self.num_dupl > 1:
@@ -420,6 +424,9 @@ def set_variational_forms(self):
         self.deltaW_kin, self.deltaW_kin_old, self.deltaW_kin_mid = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0)
         self.deltaW_int, self.deltaW_int_old, self.deltaW_int_mid = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0)
         self.deltaW_p,   self.deltaW_p_old,   self.deltaW_p_mid   = [], [], []
+        # element level Reynolds number components
+        self.Re_c,      self.Re_c_old,        self.Re_c_mid       = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0)
+        self.Re_ktilde, self.Re_ktilde_old,   self.Re_ktilde_mid  = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0)
 
         for n, M in enumerate(self.domain_ids):
 
@@ -454,6 +461,15 @@ def set_variational_forms(self):
             self.deltaW_p_old.append( self.vf.deltaW_int_pres(self.v_old, self.var_p_[j], self.io.dx(M), F=self.alevar['Fale_old']) )
             self.deltaW_p_mid.append( self.vf.deltaW_int_pres(self.vel_mid, self.var_p_[j], self.io.dx(M), F=self.alevar['Fale_mid']) )
 
+            # element level Reynolds number components - not used so far! Need to assemble a cell-based vector in order to evaluate these...
+            self.Re_c      += self.vf.re_c(self.rho[n], self.v, self.io.dx(M), w=self.alevar['w'], F=self.alevar['Fale'])
+            self.Re_c_old  += self.vf.re_c(self.rho[n], self.v_old, self.io.dx(M), w=self.alevar['w_old'], F=self.alevar['Fale_old'])
+            self.Re_c_mid  += self.vf.re_c(self.rho[n], self.vel_mid, self.io.dx(M), w=self.alevar['w_mid'], F=self.alevar['Fale_mid'])
+
+            self.Re_ktilde     += self.vf.re_ktilde(self.rho[n], self.v, self.io.dx(M), w=self.alevar['w'], F=self.alevar['Fale'])
+            self.Re_ktilde_old += self.vf.re_ktilde(self.rho[n], self.v_old, self.io.dx(M), w=self.alevar['w_old'], F=self.alevar['Fale_old'])
+            self.Re_ktilde_mid += self.vf.re_ktilde(self.rho[n], self.vel_mid, self.io.dx(M), w=self.alevar['w_mid'], F=self.alevar['Fale_mid'])
+
         # external virtual power (from Neumann or Robin boundary conditions, body forces, ...)
         w_neumann, w_body, w_robin, w_stabneumann, w_robin_valve, w_membrane = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0)
         w_neumann_old, w_body_old, w_robin_old, w_stabneumann_old, w_robin_valve_old, w_membrane_old = ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0), ufl.as_ufl(0)
@@ -553,6 +569,13 @@ def set_variational_forms(self):
             assert(self.order_vel==self.order_pres)
 
             vscale = self.stabilization['vscale']
+            try: vscale_vel_dep = self.stabilization['vscale_vel_dep']
+            except: vscale_vel_dep = False
+
+            if vscale_vel_dep:
+                vscale_max = ufl.max_value(ufl.sqrt(ufl.dot(self.v_old,self.v_old)), vscale)
+            else:
+                vscale_max = vscale
 
             h = self.io.hd0 # cell diameter (could also use max edge length self.io.emax0, but seems to yield similar/same results)
 
@@ -567,9 +590,9 @@ def set_variational_forms(self):
                     if self.num_dupl==1: j=0
                     else: j=n
 
-                    tau_supg = h / vscale
-                    tau_pspg = h / vscale
-                    tau_lsic = h * vscale
+                    tau_supg = h / vscale_max
+                    tau_pspg = h / vscale_max
+                    tau_lsic = h * vscale_max
 
                     # strong momentum residuals
                     if self.fluid_governing_type=='navierstokes_transient':
@@ -593,13 +616,13 @@ def set_variational_forms(self):
 
                     # SUPG (streamline-upwind Petrov-Galerkin) for Navier-Stokes
                     if self.fluid_governing_type=='navierstokes_transient' or self.fluid_governing_type=='navierstokes_steady':
-                        self.deltaW_int     += self.vf.stab_supg(self.acc, self.v, self.p_[j], residual_v_strong, tau_supg, self.rho[n], self.io.dx(M), w=self.alevar['w'], F=self.alevar['Fale'], symmetric=symm)
-                        self.deltaW_int_old += self.vf.stab_supg(self.a_old, self.v_old, self.p_old_[j], residual_v_strong_old, tau_supg, self.rho[n], self.io.dx(M), w=self.alevar['w_old'], F=self.alevar['Fale_old'], symmetric=symm)
-                        self.deltaW_int_mid += self.vf.stab_supg(self.acc_mid, self.vel_mid, self.pf_mid_[j], residual_v_strong_mid, tau_supg, self.rho[n], self.io.dx(M), w=self.alevar['w_mid'], F=self.alevar['Fale_mid'], symmetric=symm)
+                        self.deltaW_int     += self.vf.stab_supg(self.v, residual_v_strong, tau_supg, self.rho[n], self.io.dx(M), w=self.alevar['w'], F=self.alevar['Fale'], symmetric=symm)
+                        self.deltaW_int_old += self.vf.stab_supg(self.v_old, residual_v_strong_old, tau_supg, self.rho[n], self.io.dx(M), w=self.alevar['w_old'], F=self.alevar['Fale_old'], symmetric=symm)
+                        self.deltaW_int_mid += self.vf.stab_supg(self.vel_mid, residual_v_strong_mid, tau_supg, self.rho[n], self.io.dx(M), w=self.alevar['w_mid'], F=self.alevar['Fale_mid'], symmetric=symm)
                     # PSPG (pressure-stabilizing Petrov-Galerkin) for Navier-Stokes and Stokes
-                    self.deltaW_p[n]     += self.vf.stab_pspg(self.acc, self.v, self.p_[j], self.var_p_[j], residual_v_strong, tau_pspg, self.rho[n], self.io.dx(M), F=self.alevar['Fale'])
-                    self.deltaW_p_old[n] += self.vf.stab_pspg(self.a_old, self.v_old, self.p_old_[j], self.var_p_[j], residual_v_strong_old, tau_pspg, self.rho[n], self.io.dx(M), F=self.alevar['Fale_old'])
-                    self.deltaW_p_mid[n] += self.vf.stab_pspg(self.acc_mid, self.vel_mid, self.pf_mid_[j], self.var_p_[j], residual_v_strong_mid, tau_pspg, self.rho[n], self.io.dx(M), F=self.alevar['Fale_mid'])
+                    self.deltaW_p[n]     += self.vf.stab_pspg(self.var_p_[j], residual_v_strong, tau_pspg, self.rho[n], self.io.dx(M), F=self.alevar['Fale'])
+                    self.deltaW_p_old[n] += self.vf.stab_pspg(self.var_p_[j], residual_v_strong_old, tau_pspg, self.rho[n], self.io.dx(M), F=self.alevar['Fale_old'])
+                    self.deltaW_p_mid[n] += self.vf.stab_pspg(self.var_p_[j], residual_v_strong_mid, tau_pspg, self.rho[n], self.io.dx(M), F=self.alevar['Fale_mid'])
                     # LSIC (least-squares on incompressibility constraint) for Navier-Stokes and Stokes
                     self.deltaW_int     += self.vf.stab_lsic(self.v, tau_lsic, self.rho[n], self.io.dx(M), F=self.alevar['Fale'])
                     self.deltaW_int_old += self.vf.stab_lsic(self.v_old, tau_lsic, self.rho[n], self.io.dx(M), F=self.alevar['Fale_old'])
@@ -620,9 +643,9 @@ def set_variational_forms(self):
                     dscales = self.stabilization['dscales']
 
                     # yields same result as above for navierstokes_steady
-                    delta1 = dscales[0] * self.rho[n] * h / vscale
-                    delta2 = dscales[1] * self.rho[n] * h * vscale
-                    delta3 = dscales[2] * h / vscale
+                    delta1 = dscales[0] * self.rho[n] * h / vscale_max
+                    delta2 = dscales[1] * self.rho[n] * h * vscale_max
+                    delta3 = dscales[2] * h / vscale_max
 
                     self.deltaW_int     += self.vf.stab_v(delta1, delta2, delta3, self.v, self.p_[j], self.io.dx(M), w=self.alevar['w'], F=self.alevar['Fale'], symmetric=symm)
                     self.deltaW_int_old += self.vf.stab_v(delta1, delta2, delta3, self.v_old, self.p_old_[j], self.io.dx(M), w=self.alevar['w_old'], F=self.alevar['Fale_old'], symmetric=symm)
@@ -632,6 +655,7 @@ def set_variational_forms(self):
                     self.deltaW_p_old[n] += self.vf.stab_p(delta1, delta3, self.v_old, self.p_old_[j], self.var_p_[j], self.rho[n], self.io.dx(M), w=self.alevar['w_old'], F=self.alevar['Fale_old'])
                     self.deltaW_p_mid[n] += self.vf.stab_p(delta1, delta3, self.vel_mid, self.pf_mid_[j], self.var_p_[j], self.rho[n], self.io.dx(M), w=self.alevar['w_mid'], F=self.alevar['Fale_mid'])
 
+
         ### full weakforms
 
         # kinetic plus internal minus external virtual power

src/ambit_fe/fluid/fluid_variationalform.py

@@ -61,28 +61,44 @@ def deltaW_int_pres(self, v, var_p, ddomain, w=None, F=None):
         return ufl.div(v)*var_p*ddomain
 
     def res_v_strong_navierstokes_transient(self, a, v, rho, sig, w=None, F=None):
+
+        return self.f_inert_strong_navierstokes_transient(a, v, rho, w=w, F=F) - self.f_visc_strong(sig, F=F)
+
+    def res_v_strong_navierstokes_steady(self, v, rho, sig, w=None, F=None):
+
+        return self.f_inert_strong_navierstokes_steady(v, rho, w=w, F=F) - self.f_visc_strong(sig, F=F)
+
+    def res_v_strong_stokes_transient(self, a, rho, sig, w=None, F=None):
+
+        return self.f_inert_strong_stokes_transient(a, rho) - self.f_visc_strong(sig, F=F)
+
+    def res_v_strong_stokes_steady(self, rho, sig, F=None):
+
+        return -self.f_visc_strong(sig, F=F)
+
+    def f_inert_strong_navierstokes_transient(self, a, v, rho, w=None, F=None):
         if self.formulation=='nonconservative':
-            return rho*(a + ufl.grad(v) * v) - ufl.div(sig)
+            return rho*(a + ufl.grad(v) * v)
         elif self.formulation=='conservative':
-            return rho*(a + ufl.div(ufl.outer(v,v))) - ufl.div(sig)
+            return rho*(a + ufl.div(ufl.outer(v,v)))
         else:
             raise ValueError("Unknown fluid formulation!")
 
-    def res_v_strong_navierstokes_steady(self, v, rho, sig, w=None, F=None):
+    def f_inert_strong_navierstokes_steady(self, v, rho, w=None, F=None):
         if self.formulation=='nonconservative':
-            return rho*(ufl.grad(v) * v) - ufl.div(sig)
+            return rho*(ufl.grad(v) * v)
         elif self.formulation=='conservative':
-            return rho*(ufl.div(ufl.outer(v,v))) - ufl.div(sig)
+            return rho*(ufl.div(ufl.outer(v,v)))
         else:
             raise ValueError("Unknown fluid formulation!")
 
-    def res_v_strong_stokes_transient(self, a, v, rho, sig, w=None, F=None):
+    def f_inert_strong_stokes_transient(self, a, rho):
 
-        return rho*a - ufl.div(sig)
+        return rho*a
 
-    def res_v_strong_stokes_steady(self, rho, sig, F=None):
+    def f_visc_strong(self, sig, F=None):
 
-        return -ufl.div(sig)
+        return ufl.div(sig)
 
     def res_p_strong(self, v, F=None):
 
@@ -103,14 +119,14 @@ def deltaW_ext_robin_valve(self, v, beta, dboundary, fcts='+', w=None, F=None):
 
 
     ### SUPG/PSPG stabilization - cf. Tezduyar and Osawa (2000), "Finite element stabilization parameters computed from element matrices and vectors"
-    def stab_supg(self, a, v, p, res_v_strong, tau_supg, rho, ddomain, w=None, F=None, symmetric=False):
+    def stab_supg(self, v, res_v_strong, tau_supg, rho, ddomain, w=None, F=None, symmetric=False):
 
         if symmetric: # modification to make the effective stress symmetric
             return (1./rho) * ufl.dot(tau_supg*rho*ufl.sym(ufl.grad(self.var_v))*v, res_v_strong) * ddomain
         else:
             return (1./rho) * ufl.dot(tau_supg*rho*ufl.grad(self.var_v)*v, res_v_strong) * ddomain
 
-    def stab_pspg(self, a, v, p, var_p, res_v_strong, tau_pspg, rho, ddomain, F=None):
+    def stab_pspg(self, var_p, res_v_strong, tau_pspg, rho, ddomain, F=None):
 
         return (1./rho) * ufl.dot(tau_pspg*ufl.grad(var_p), res_v_strong) * ddomain
 
@@ -136,6 +152,16 @@ def stab_p(self, delta1, delta3, v, p, var_p, rho, ddomain, w=None, F=None):
         return (1./rho) * ( delta1 * ufl.dot(ufl.grad(v)*v, ufl.grad(var_p)) + \
                             delta3 * ufl.dot(ufl.grad(p), ufl.grad(var_p)) ) * ddomain
 
+    # components of element-level Reynolds number - cf. Tezduyar and Osawa (2000)
+    def re_c(self, rho, v, ddomain, w=None, F=None):
+
+        return rho * ufl.dot(ufl.grad(v)*v, self.var_v) * ddomain
+
+    def re_ktilde(self, rho, v, ddomain, w=None, F=None):
+
+        return rho * ufl.dot(ufl.grad(v)*v, ufl.grad(self.var_v)*v) * ddomain
+
+
     ### Flux coupling conditions
 
     # flux
@@ -242,33 +268,49 @@ def deltaW_int_robin_cur(self, v, vD, beta, dboundary, F=None, fcts=None):
             return (beta*ufl.dot((v-vD), self.var_v) * J*ufl.sqrt(ufl.dot(self.n0, (ufl.inv(F)*ufl.inv(F).T)*self.n0)))(fcts)*dboundary
 
     def res_v_strong_navierstokes_transient(self, a, v, rho, sig, w=None, F=None):
-        i, j, k = ufl.indices(3)
+
+        return self.f_inert_navierstokes_transient(a, v, rho, w=w, F=F) - self.f_visc_strong(sig, F=F)
+
+    def res_v_strong_navierstokes_steady(self, v, rho, sig, w=None, F=None):
+
+        return self.f_inert_navierstokes_steady(v, rho, w=w, F=F) - self.f_visc_strong(sig, F=F)
+
+    def res_v_strong_stokes_transient(self, a, v, rho, sig, w=None, F=None):
+
+        return self.f_inert_stokes_transient(a, v, rho, w=w, F=F) - self.f_visc_strong(sig, F=F)
+
+    def res_v_strong_stokes_steady(self, rho, sig, F=None):
+
+        return -self.f_visc_strong(sig, F=F)
+
+    def f_inert_navierstokes_transient(self, a, v, rho, sig, w=None, F=None):
         if self.formulation=='nonconservative':
-            return rho*(a + ufl.grad(v)*ufl.inv(F) * (v-w)) - ufl.as_vector(ufl.grad(sig)[i,j,k]*ufl.inv(F).T[j,k], i)
+            return rho*(a + ufl.grad(v)*ufl.inv(F) * (v-w))
         elif self.formulation=='conservative':
-            return rho*(a + ufl.as_vector(ufl.grad(ufl.outer(v,v-w))[i,j,k]*ufl.inv(F).T[j,k], i)) - ufl.as_vector(ufl.grad(sig)[l,m,n]*ufl.inv(F).T[j,k], i)
+            i, j, k = ufl.indices(3)
+            return rho*(a + ufl.as_vector(ufl.grad(ufl.outer(v,v-w))[i,j,k]*ufl.inv(F).T[j,k], i))
         else:
             raise ValueError("Unknown fluid formulation!")
 
-    def res_v_strong_navierstokes_steady(self, v, rho, sig, w=None, F=None):
-        i, j, k = ufl.indices(3)
+    def f_inert_navierstokes_steady(self, v, rho, sig, w=None, F=None):
         if self.formulation=='nonconservative':
-            return rho*(ufl.grad(v)*ufl.inv(F) * (v-w)) - ufl.as_vector(ufl.grad(sig)[i,j,k]*ufl.inv(F).T[j,k], i)
+            return rho*(ufl.grad(v)*ufl.inv(F) * (v-w))
         elif self.formulation=='conservative':
-            return rho*(ufl.as_vector(ufl.grad(ufl.outer(v,v-w))[i,j,k]*ufl.inv(F).T[j,k], i)) - ufl.as_vector(ufl.grad(sig)[l,m,n]*ufl.inv(F).T[j,k], i)
+            i, j, k = ufl.indices(3)
+            return rho*(ufl.as_vector(ufl.grad(ufl.outer(v,v-w))[i,j,k]*ufl.inv(F).T[j,k], i))
         else:
             raise ValueError("Unknown fluid formulation!")
 
-    def res_v_strong_stokes_transient(self, a, v, rho, sig, w=None, F=None):
-        i, j, k = ufl.indices(3)
+    def f_inert_stokes_transient(self, a, v, rho, w=None, F=None):
         if self.formulation=='nonconservative':
-            return rho*(a + ufl.grad(v)*ufl.inv(F) * (-w)) - ufl.as_vector(ufl.grad(sig)[i,j,k]*ufl.inv(F).T[j,k], i)
+            return rho*(a + ufl.grad(v)*ufl.inv(F) * (-w))
         elif self.formulation=='conservative':
-            return rho*(a + ufl.as_vector(ufl.grad(ufl.outer(v,-w))[i,j,k]*ufl.inv(F).T[j,k], i)) - ufl.as_vector(ufl.grad(sig)[i,j,k]*ufl.inv(F).T[j,k], i)
+            i, j, k = ufl.indices(3)
+            return rho*(a + ufl.as_vector(ufl.grad(ufl.outer(v,-w))[i,j,k]*ufl.inv(F).T[j,k], i))
         else:
             raise ValueError("Unknown fluid formulation!")
 
-    def res_v_strong_stokes_steady(self, rho, sig, F=None):
+    def f_visc_strong(self, sig, F=None):
         i, j, k = ufl.indices(3)
         return -ufl.as_vector(ufl.grad(sig)[i,j,k]*ufl.inv(F).T[j,k], i)
 
@@ -291,14 +333,14 @@ def deltaW_ext_robin_valve(self, v, beta, dboundary, fcts='+', w=None, F=None):
 
 
     ### SUPG/PSPG stabilization
-    def stab_supg(self, a, v, p, res_v_strong, tau_supg, rho, ddomain, w=None, F=None, symmetric=False):
+    def stab_supg(self, v, res_v_strong, tau_supg, rho, ddomain, w=None, F=None, symmetric=False):
         J = ufl.det(F)
         if symmetric: # modification to make the effective stress symmetric
             return (1./rho) * ufl.dot(tau_supg*rho*ufl.sym(ufl.grad(self.var_v)*ufl.inv(F))*v, res_v_strong) * J*ddomain
         else:
             return (1./rho) * ufl.dot(tau_supg*rho*ufl.grad(self.var_v)*ufl.inv(F)*v, res_v_strong) * J*ddomain
 
-    def stab_pspg(self, a, v, p, var_p, res_v_strong, tau_pspg, rho, ddomain, F=None):
+    def stab_pspg(self, var_p, res_v_strong, tau_pspg, rho, ddomain, F=None):
         J = ufl.det(F)
         return (1./rho) * ufl.dot(tau_pspg*ufl.inv(F).T*ufl.grad(var_p), res_v_strong) * J*ddomain
 
@@ -322,6 +364,15 @@ def stab_p(self, delta1, delta3, v, p, var_p, rho, ddomain, w=None, F=None):
         return (1./rho) * ( delta1 * ufl.dot(ufl.grad(v)*ufl.inv(F)*(v-w), ufl.inv(F).T*ufl.grad(var_p)) + \
                             delta3 * ufl.dot(ufl.inv(F).T*ufl.grad(p), ufl.inv(F).T*ufl.grad(var_p)) ) * J*ddomain
 
+    # components of element-level Reynolds number
+    def re_c(self, rho, v, ddomain, w=None, F=None):
+        J = ufl.det(F)
+        return rho * ufl.dot(ufl.grad(v)*ufl.inv(F)*(v-w), self.var_v) * J*ddomain
+
+    def re_ktilde(self, rho, v, ddomain, w=None, F=None):
+        J = ufl.det(F)
+        return rho * ufl.dot(ufl.grad(v)*ufl.inv(F)*(v-w), ufl.grad(self.var_v)*ufl.inv(F)*v) * J*ddomain
+
 
     ### Flux coupling conditions
 

src/ambit_fe/ioroutines.py

@@ -900,6 +900,10 @@ def write_output(self, pb, writemesh=False, N=1, t=0):
                         pw_out = fem.Function(pb.V_out_scalar, name=pw.name)
                         pw_out.interpolate(pw)
                         self.resultsfiles[res].write_function(pw_out, indicator)
+                    elif res=='reynolds':
+                        re_out = fem.Function(pb.V_out_vector, name=pb.re.name)
+                        re_out.interpolate(pb.re)
+                        self.resultsfiles[res].write_function(re_out, indicator)
                     elif res=='counters':
                         # iteration counters, written by base class
                         pass

src/ambit_fe/solver/preconditioner.py

@@ -521,7 +521,7 @@ def setUp(self, pc):
         self.Et.copy(result=self.Tmod)
         self.Tmod.axpy(-1., self.B_Adinv_Dt)
 
-        # --- Wmod = R - D diag(A)^{-1} Dt - E diag(Smod)^{-1} Tmod + D diag(A)^{-1} Bt diag(Smod)^{-1} Tmod
+        # --- Wmod = R - D diag(A)^{-1} Dt - Umod diag(Smod)^{-1} Tmod
 
         if self.schur_block_scaling[1]['type']=='diag':
             self.Smod.getDiagonal(result=self.smoddinv_vec)