Help with Building Block Matrix

General Discussion
1

Hi everyone,

I’m currently working on finding a solution to Stokes flow within a square domain. Here’s a brief overview of the problem:

Domain Shape: Square-like, with the top and bottom boundaries representing no-penetration and no-slip conditions, respectively.

Boundary Conditions: The left and right boundaries are set as periodic inlets/outlets.
Constraint: I need to impose a constraint where the fluid flux along the left/right boundaries is prescribed by a user-defined number (let’s arbitrarily set it to 1).

I’ve been trying to tackle this by using a block matrix solve, but I’m encountering an error during the construction of my A matrix (which represents the Stokes equation). Despite spending considerable time on it, I haven’t been able to make much headway.

I would greatly appreciate any help, whether it’s providing alternative strategies or suggesting fixes to my current approach. Thank you in advance for your assistance!

-Fluidman

StokesFlux.edp (2.2 KB)

2

This looks like a problem with the array syntax. I defined an fespace like this
and changed your varf to get the matrix but didn’t have time to fix the rhs
see if you can get this to work,

fespace Xh(Th, [P2,P2,P2], periodic=[[Perl, y],[Perr, y]]);
3

Hi Mike,

Thanks for your help earlier.

After incorporating your suggestions, I continued to refine the code, and now it’s executing properly! However, I believe there’s a conceptual gap in my understanding because the solution I’m obtaining is only the size of a single one of my unknowns. In other words, when I compute x = (AA^{-1})b, where AA is my block matrix, the resulting x vector is only the size of a single one of my unknowns (plus 1 due to the Lagrange multiplier). I suspect I’m overlooking a key aspect of my problem. Below is the updated code, including print commands to display the sizes of the various matrices and vectors involved.
-Fluid Man
StokesFlux.edp (2.8 KB)

4

I always have a hard time with the vector or array notation since I want
to think about how the dof’s are organized in memory and in the varf parameter
list you can leave them out and its splits the list into trial and test functions as scalars’
However that can be awkward too.

IIRC referencing the first component is like referncing the whole thing or something
similar- hopefully someone else can comment :slight_smile:

You can see this more easily by removing the periodic BC and change to P1
elements ( often that is defined as a macro PK P1 altough they need not all be the same ).

If you do this, you can see the number of verticies is 1/3 the array size. I haven’t
written any FF code in a while and it took a while :slight_smile:

5

Dear Dylan,
Since you defined
fespace Vh(Th, [P2,P2,P2], periodic=[[Perl, y],[Perr, y]]);
Vh [ux, uy, p];
you have to understand [ux,uy,p] as a full vector of unknowns. The thing that is difficult to understand is that FreeFem calls the list of its components in the associated base as ux[] so here it means the whole set of ux,uy,p.
Other equivalent names of the same vector are uy[] or p[].
It you want to separate the components you have to set

fespace Gh(Th, P2, periodic=[[Perl, y],[Perr, y]]);
Gh ux1,uy1,p1;
ux1=ux;
uy1=uy;
p1=p;

Then ux1 is a separate FE function with size the dimension of Gh, and ux1[] is the vector of the components of ux. The same statement is valid for uy1 and p1.
Here is your code with little corrections.
StokesFlux.edp (2.7 KB)

Francois.