How to assemble the right-hand side vector?

[liqihao2000]

Dear Mikael,
I'm constructing a decoupling algorithm for a Cahn-Hilliard Navier-Stokes system.

$$ \phi_t + \nabla\cdot({\bf u}\phi) = M \Delta \mu,$$

$$ \mu = \lambda(-\Delta \phi) + f(\phi)), $$

$$ {\bf u}_t + ({\bf u}\cdot\nabla){\bf u} -\nu\Delta{\bf u} + \nabla p + \phi\nabla \mu =0, $$

$$ \nabla\cdot {\bf u} =0. $$

How to assemble the right-hand side vector
$$-\int_\Omega \psi \nabla\cdot({\bf u}^n\phi^n) ~dx ,$$
where the ${\bf u}^n$ is known which is the initial velocity vector in the Navier-Stokes equation, $\phi^n$ and $\mu^n$ are known in Cahn-Hilliard equation, $\psi$ is the test function for the Cahn-Hilliard equation.

The minimal working example (MWE) is as follows:

import matplotlib.pyplot as plt
import sympy as sp
from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, FunctionSpace, TensorSpace, \
    TensorProductSpace, VectorSpace, CompositeSpace, BlockMatrix, Array, comm, Dx, dot, outer, project

N = (64, 64)

# basis function for velocity components in x and y directions: P_{N}
bcx = {'left': {'N': 0}, 'right': {'N': 0}}  # BCs of phi
bcy = {'left': {'N': 0}, 'right': {'N': 0}}  # BCs of phi
K0 = FunctionSpace(N[0], family='Legendre', bc=bcx, domain=(-1, 1))  # function space for phi
K1 = FunctionSpace(N[1], family='Legendre', bc=bcy, domain=(-1, 1))  # function space for phi
K0X = FunctionSpace(N[0], family='Legendre', bc=(0, 0), domain=(-1, 1))  # for velocity
K0Y = FunctionSpace(N[1], family='Legendre', bc=(0, 0), domain=(-1, 1))  # for velocity

# basis function for pressure: P_{N-2}
PX = FunctionSpace(N[0], 'Legendre', quad='GL', domain=(-1, 1))
PY = FunctionSpace(N[1], 'Legendre', quad='GL', domain=(-1, 1))
PX.slice = lambda: slice(0, N[0] - 2)
PY.slice = lambda: slice(0, N[1] - 2)

# Define a multi-dimensional tensor product basis
Vs = TensorProductSpace(comm, (K0X, K0Y))  # velocity in x-direction
V0 = Vs.get_homogeneous()  # velocity in y-direction

phi_function_space = TensorProductSpace(comm, (K0, K1), axes=(0, 1))
CH_function_space = CompositeSpace([phi_function_space, phi_function_space])  # mixed function space for CH
grad_phi_function_space = VectorSpace([phi_function_space, phi_function_space])
pressure_function_space = TensorProductSpace(comm, (PX, PY), modify_spaces_inplace=True)

# Create vector space for velocity
velocity_function_space = VectorSpace([Vs, V0])  # u1
velocity2_function_space = VectorSpace([V0, V0])  # u2
stress_function_space = TensorSpace([velocity_function_space, velocity_function_space])  # cauchy stress tensor

# Parameter
lambd = 0.1
epsilon = 0.1

x, y, t = sp.symbols("x,y,t", real=True)

ex_phi = sp.cos(2 * sp.pi * y) * sp.cos(2 * sp.pi * x) * sp.sin(t)
ex_mu = lambd * (- (ex_phi.diff(x, 2) + ex_phi.diff(y, 2)) + (ex_phi ** 3 - ex_phi) / (epsilon ** 2))
ex_u = (sp.sin(sp.pi * x) ** 2) * sp.sin(2 * sp.pi * y) * sp.sin(t)
ex_v = -sp.sin(2 * sp.pi * x) * (sp.sin(sp.pi * y) ** 2) * sp.sin(t)
ex_p = sp.cos(sp.pi * x) * sp.cos(sp.pi * y) * sp.sin(t)

# Create test and trial spaces for phi, mu, velocity and pressure
phi_trial, mu_trial = TrialFunction(CH_function_space)
psi_test, omega_test = TestFunction(CH_function_space)
u_trial, v_test = TrialFunction(velocity_function_space), TestFunction(velocity_function_space)

time = 0.1
phi_prev = Array(phi_function_space, buffer=ex_phi.subs(t, time))
mu_prev = Array(phi_function_space, buffer=ex_mu.subs(t, time))
u_prev = Array(velocity_function_space, buffer=(ex_u.subs(t, time), ex_v.subs(t, time)))
p_prev = Array(pressure_function_space, buffer=ex_p.subs(t, time))

rhs_sub3 = Function(CH_function_space)
rhs1_sub3, rhs2_sub3 = rhs_sub3

rhs_velocity_u2 = Function(velocity_function_space)

# Compute right and side
uphi = u_prev * phi_prev
uphi_hat = uphi.forward()

# Assemble right hand side vector

How to assemble the right-hand side vector
$$-\int_\Omega \psi \nabla\cdot({\bf u}^n\phi^n) ~dx ,$$

# Question 1 ? 
div_uphi_hat = project(div(uphi_hat), Vs)
rhs1_sub3 = inner(psi_test, -div_uphi_hat.backward(), output_array=rhs1_sub3)

[mikaem]

Hi @liqihao2000
I'm sorry, I don't know how I could have missed your question. Have you sorted this out yet?

[liqihao2000]

Hi Mikael, it doesn't matter. I don't know how to get divergence, so I haven't solved this problem yet.

[mikaem]

Hi
You quite correctly get the divergence by projection, but you probably need to use a space with no associated boundary conditions

div_uphi_hat = project(div(uphi_hat), Vs.get_orthogonal())

Spaces with boundary conditions are generally only for solving PDEs and should normally not be used with projections.

Did you see this demo? https://github.com/spectralDNS/shenfun/blob/master/demo/NavierStokesDrivenCavity.py
There are some examples of computing similar terms there.