I'm trying to implement the paper "Conformation Constraints for Efficient Viscoelastic Fluid Simulation": https://www.gmrv.es/Publications/2017/BGAO17/SA2017_BGAO.pdf

In implementing eqn. 6, they solve the following linear system for \(Q\):

$$
\frac{Q - Q_0}{\Delta t} = Q \nabla \vec{u} + (\nabla \vec{u})^T Q - \frac{1}{\tau} (Q - I)
$$

\(Q_0, \nabla \vec{u}, \Delta t\) and \(\tau\) are all knowns. \(Q\) is a symmetric 3x3 matrix.

In the paper, they describe reformulating \(Q\) as a 6-component vector: $$q = [Q_{xx},Q_{yy},Q_{zz},Q_{xy},Q_{xz},Q_{yz}]^T$$

With \(I\) becoming the vector: $$I = \bar{q} = [1,1,1,0,0,0]^T$$

To be able to solve this in practice I think I need to get it in the form \(Aq = b\).

I'd love some help figuring out how to reformulate the above equation into a matrix-vector product that I can use a linear solver on, as I can't seem to just naïvely rearrange terms.

You should probably ask this kind of questions in a numerics or fluid dynamics forum. Anyways, I don't have enough time to read the paper but maybe this helps a little bit:

Let's treat Q as a 3x3 matrix. So all additive terms have to be 3x3 matrices too. Now remove the parenthesis:

1/dt Q - 1/dt Q0 =  Q grad(u) + grad(u)TQ - 1/tau Q - 1/tau I

separate dependent and independent terms:

1/tau I - 1/dt Q0 =  Q grad(u) + grad(u)TQ - 1/tau Q - 1/dt Q

Substitution:

R = 1/tau I - 1/dt Q0

R =  Q grad(u) + grad(u)TQ - 1/tau Q - 1/dt Q

R =  Q grad(u) + grad(u)TQ - (1/tau + 1/dt) Q

Now with

| q11, q12, q13 |    | q11, q12, q13 |

| q21, q22, q23 | = | q12, q22, q23 | = Q

| q31, q32, q33 |    | q13, q23, q33 |

and

| du1_dx, du2_dx, du3_dx |

| du1_dy, du2_dy, du3_dy | = grad(u)   --- (Not sure if it is not grad(u)T instead)

| du1_dz, du2_dz, du3_dz |

you should explicitly calculate the result of

Q grad(u) + grad(u)TQ

By doing so, you can probably figure out, how the terms need to be arranged in a matrix A to rewrite the term as

A q

where q is now a 6x1 vector.

Transforming all other terms to the vector notation should give you an equation like this:

r =  A q - (1/tau + 1/dt) q

Now "smuggle in" an Identity matrix:

r =  A q - (1/tau + 1/dt) I q

Then you get something like this:

r =  (A - (1/tau + 1/dt) I) q

This can now be solved for q with any kind of linear solver. No guarantees that everything is correct here. For the substep when you transform Q grad(u) + grad(u)TQ it might be helpful to google something like "stiffness matrix voigt notation" because it is a quite similar problem. If you still have problems, feel free to ask, even though I don't have too much time at the moment and it's not exactly a topic where I have a lot of experience.

Greetings