Oldroyd-B model freefem++ code bugs and issues: Assertion fail

ozil101 · February 24, 2025, 11:19am

Hi all,

I built up a code for 2D bearing lubrication problem, I adjusted all the parameters & conditions into the standard ones of this benchmark problem, however I always get this issue:

{ 139 : } sizestack + 1024 =20004 ( 18980 )

– mesh: Nb of Triangles = 4382, Nb of Vertices 2316
current line = 50
Assertion fail : (!half)
line :1263, in file …/femlib/HashMatrix.cpp
Assertion fail : (!half)
line :1263, in file …/femlib/HashMatrix.cpp
err code 6 , mpirank 0}

what could be the reason and the solution?

best,
Hussein.

frederichecht · February 24, 2025, 12:27pm

Without the code it is difficult,It you say the building matrix is full but you store just half.
So You do strange thing.

ozil101 · February 24, 2025, 3:16pm

Thank you Prof for your reply, this is the code. I tried to implement the conditions within this paper [https://onlinelibrary.wiley.com/doi/10.1002/zamm.202300095
Bmmpff_bearing_lubrication.edp (4.8 KB)
].

frederichecht · February 24, 2025, 3:40pm

I have run you script , I a see a first mistake the boundary condition on Oldroyd… are take on table 4 but the label 4 do not existe => full Neumann Boundary condition.

ozil101 · February 27, 2025, 12:02pm

Thank you the code works now very well. However, I COULD NOT manage to apply some FEM stabilization techniques (i.e. SUPG, PGP, DEVSS) on this benchmark I always got an operator error especially when applying DEVSS (using least-square scheme by adding auxiliary velocity gradient variables: (G11, G12, G21, G22)) scheme like this:
{
{ 105 : solve ComputeG(G11, G12, G21, G22) =
106 : int2d(Th)(G11 * w - dx(u) * w) error operator - <10LinearCombI7MGauche4C_F0E>,
List of choices
( : , )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <14Matrice_CreuseISt7complexIdEE>, <14Matrice_CreuseISt7complexIdEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <14Matrice_CreuseISt7complexIdEE>, <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE>, <14Matrice_CreuseISt7complexIdEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE>, <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <14Matrice_CreuseIdE>, <14Matrice_CreuseIdE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <14Matrice_CreuseIdE>, <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE>, <14Matrice_CreuseIdE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE>, <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> )
( <10FormLinear> : <10FormLinear> )
( <12FormBilinear> : <12FormBilinear> )
( <6C_args> : <6C_args>, <10FormLinear> )
( <6C_args> : <6C_args>, <12FormBilinear> )
( <10LinearCombI7MGauche4C_F0E> : <10LinearCombI7MGauche4C_F0E>, <10LinearCombI7MGauche4C_F0E> )
( <10LinearCombI6MDroit4C_F0E> : <10LinearCombI6MDroit4C_F0E>, <10LinearCombI6MDroit4C_F0E> )
( <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> : <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E>, <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> )
( <10LinearCombI7MGauche4C_F0E> : <10LinearCombI7MGauche4C_F0E> )
( <10LinearCombI6MDroit4C_F0E> : <10LinearCombI6MDroit4C_F0E> )
( <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> : <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> )
( <4C_F0> : <12TransE_Array>, <12TransE_Array> )
( <4C_F0> : <7E_Array>, <12TransE_Array> )
( <4C_F0> : <7E_Array>, <7E_Array> )
( <4C_F0> : <12TransE_Array>, <7E_Array> )
( <12Add_Mulc_KN_ISt7complexIdEE> : <8Mulc_KN_ISt7complexIdEE>, <3KN_ISt7complexIdEE> )
( <12Add_Mulc_KN_ISt7complexIdEE> : <3KN_ISt7complexIdEE>, <8Mulc_KN_ISt7complexIdEE> )
( <12Add_Mulc_KN_ISt7complexIdEE> : <8Mulc_KN_ISt7complexIdEE>, <8Mulc_KN_ISt7complexIdEE> )
( <8Mulc_KN_ISt7complexIdEE> : <3KN_ISt7complexIdEE> )
( <7Sub_KN_ISt7complexIdEE> : <3KN_ISt7complexIdEE>, <3KN_ISt7complexIdEE> )
( <12Add_Mulc_KN_IdE> : <8Mulc_KN_IdE>, <3KN_IdE> )
( <12Add_Mulc_KN_IdE> : <3KN_IdE>, <8Mulc_KN_IdE> )
( <12Add_Mulc_KN_IdE> : <8Mulc_KN_IdE>, <8Mulc_KN_IdE> )
( <8Mulc_KN_IdE> : <3KN_IdE> )
( <7Sub_KN_IdE> : <3KN_IdE>, <3KN_IdE> )
( <12Add_Mulc_KN_IlE> : <8Mulc_KN_IlE>, <3KN_IlE> )
( <12Add_Mulc_KN_IlE> : <3KN_IlE>, <8Mulc_KN_IlE> )
( <12Add_Mulc_KN_IlE> : <8Mulc_KN_IlE>, <8Mulc_KN_IlE> )
( <8Mulc_KN_IlE> : <3KN_IlE> )
( <7Sub_KN_IlE> : <3KN_IlE>, <3KN_IlE> )
( : )
( : )
( : )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )

Error line number 106, in file …edp, before token )

current line = 106
Compile error :
line number :106, )
error Compile error :
line number :106, )
code = 1 mpirank: 0
}

frederichecht · February 27, 2025, 12:31pm

Remember, I, freefem with solve variation form with unknown function and test function, and in the example no test function so we get an error!

ozil101 · February 27, 2025, 1:10pm

Yes, indeed Prof. I also re-define my constitutive Oldroyd-B model solution as this and also defined the initial variable in my fe-space and test function (I use ‘w’ test function as the same old square-root conformation solution and defined G11old, G11, G12old, G12, G13old, G21, G14old, G22):
Full code:
{
int outer = 1;
int inner = 2;

border a(t=0, 2pi){x=cos(t); y=sin(t); label=outer;}
border b(t=0, 2pi){x=0.55+0.4cos(t); y=0.35sin(t); label=inner;}

plot(a(50) + b(30));

mesh Th = buildmesh(a(350) + b(-300));
plot(Th, fill=true, ps=“Bmmpff_Bearing_Lubrication_Domain.eps”);

fespace Vh(Th, P1);
fespace Mh(Th, P1);

real dt = 0.005, re = 500.0, Wi = 2.75;
real al=0.5, dT=dt/2.0, t=0, T=30.0, ee=0.00001, L2errortau11, L2errortau12, L2erroru;
real gamma = 0.1, mur = 0.50 / re, sigmanf = 0.5, B0 = 1.1;
real rhonf = 0.25, j = 0.1;
real munf = 0.525; // fluid dynamic viscosity value
real Omega = 0.1025; // initial fluid angular velocity variable
real lambda = 0.505; // extra stress tensor parameter
real beta=1.0/9.0;
real epsilon = 0.314;

real nv = Th.nv;
real nt = Th.nt;

Mh p, q;
Vh w, u, v, omega, tau11, tau12, tau21, tau22, G11, G12, G21, G22;
Vh uT, vT, omegaT, uTT, vTT, omegaTT;
Vh tau11old, tau12old, tau21old, tau22old, uold, vold, pold, omegaold, G11old, G12old, G13old, G14old;

u = -0.5 * Omega * sin(epsilon);
v = 0.5 * Omega * cos(epsilon);
omega = 200;

for (int m = 0; m <= T / dt; m++) {
Vh uold = u, vold = v, pold = p, omegaold = omega;
pold = p;

//Intermediate Velocity Solutions//

solve Adual4unon(uT, w) =
    int2d(Th)(rhonf * uT * w / dt)
    + int2d(Th)(munf *  (dx(uT) * dx(w) + dy(uT) * dy(w)))
    - int2d(Th)((rhonf * convect([uold, vold], -dt, uold) * w) / dt)
    + int2d(Th)(2 * rhonf * Omega * (dy(vold) * w - dx(uold) * w))
    - int2d(Th)(sigmanf * B0^2 * uT * w)
    + int2d(Th)(lambda * (dx(uold) * dx(w) + dy(uold) * dy(w)))  
    + on(outer, uT = -Omega * y)
    + on(inner, uT = 0);  // No flow inside the small circle

solve Adual4vnon(vT, w) =
    int2d(Th)(rhonf * vT * w / dt)
    + int2d(Th)(munf * (dx(vT) * dx(w) + dy(vT) * dy(w)))
    - int2d(Th)((rhonf * convect([uold, vold], -dt, vold) * w) / dt)
    + int2d(Th)(2 * rhonf * Omega * (dy(uold) * w - dx(vold) * w))
    - int2d(Th)(sigmanf * B0^2 * vT * w)
    + int2d(Th)(lambda * (dx(vold) * dx(w) + dy(vold) * dy(w)))
    + on(outer, vT = -Omega * x)
    + on(inner, vT = 0);

solve Adual4omeganon(omegaT, w) =
    int2d(Th)(omegaT * w / dt)
    - int2d(Th)(convect([uold, vold], -dt, omegaold) * w / dt)
    + on(outer, omegaT = B0)
    + on(inner, omegaT = 0);  // No microrotation inside the small circle

//Updated Velocity Solutions (u,v,omega)//

solve pb4p(p, q) =
    int2d(Th)(dx(p) * dx(q) + dy(p) * dy(q))
    + int2d(Th)(((dx(uT) + dy(vT)) * q) / dt);

solve pb4u(u, w) =
    int2d(Th)((dx(u) * dx(w) + dy(u) * dy(w)) / re + u * w / dt)
    - int2d(Th)(uT * w / dt - dx(p) * w)
    + on(outer, u = 0)
    + on(inner, u = 0);  // No flow inside the small circle

solve pb4v(v, w) =
    int2d(Th)((dx(v) * dx(w) + dy(v) * dy(w)) / re + v * w / dt)
    - int2d(Th)(vT * w / dt - dy(p) * w)
    + on(outer, v = 0)
    + on(inner, v = 0);  // No flow inside the small circle

solve pb4omega(omega, w) =
    int2d(Th)((dx(omega) * dx(w) + dy(omega) * dy(w)) / re + omega * w / dt)
    - int2d(Th)(omegaT * w / dt)
    + on(outer, omega = 0)
    + on(inner, omega = 0);  // No microrotation inside the small circle

Vh tau11old=tau11,tau12old=tau12,tau21old=tau21,tau22old=tau22;
Vh G11old=G11, G12old=G12, G13old=G21, G14old=G22;

  //**********************************solve the constitutive model equation*********************//

// define DEVSS aux. variable G
solve ComputeG(G11, G12, G21, G22) =
int2d(Th)(G11 * w - dx(u) * w)

int2d(Th)(G12 * w - dy(u) * w)
int2d(Th)(G21 * w - dx(v) * w)
int2d(Th)(G22 * w - dy(v) * w);

// include DEVSS stabilization
solve Oldroyd11 (tau11, w) =
int2d(Th)(tau11 * w / dt)

int2d(Th)(-1/dt * convect([u, v], -dt, tau11old) * w)

int2d(Th)(((1-beta) * 2 / Wi * G11 + 2 * (G11 * tau11old + G12 * tau12old) - tau11old / Wi) * w)

int1d(Th, inner)(tau11 * w)
int1d(Th, outer)(tau11 * w);

solve Oldroyd12 (tau12, w) =
int2d(Th)(tau12 * w / dt)

int2d(Th)(-1/dt * convect([u, v], -dt, tau12old) * w)

int2d(Th)(((1-beta) / Wi * (G12 + G21) + (G11 * tau12old + G12 * tau22old + G21 * tau12old + G22 * tau11old)) * w - w * tau12old / Wi)

int1d(Th, inner)(tau12 * w)
int1d(Th, outer)(tau12 * w);

tau21 = tau12;

solve Oldroyd22 (tau22, w) =
int2d(Th)(tau22 * w / dt)

int2d(Th)(-1/dt * convect([u, v], -dt, tau22old) * w)

int2d(Th)(((1-beta) * 2 / Wi * G22 + 2 * (G21 * tau12old + G22 * tau22old) - tau22old / Wi) * w)

int1d(Th, inner)(tau22 * w)
int1d(Th, outer)(tau22 * w);

t=t+dt;

// plot(u,fill=1,cmm=“t=”+t + “, min=” + u.min + “, max=” + u.max);

L2errortau11=sqrt(int2d(Th)((tau11-2.0*Wi*(1-beta)*(0.16*(1.0-2.0*y)*(1.0-2.0*y)))^2))
/sqrt(int2d(Th)((2.0*Wi*(1-beta)*(0.16*(1.0-2.0*y)*(1.0-2.0*y)))^2));
L2errortau12=sqrt(int2d(Th)((tau12-(1-beta)*(0.4*(1.0-2.0*y)))^2))/sqrt(int2d(Th)(((1-beta)*(0.4*(1.0-2.0*y)))^2));
 L2erroru=sqrt(int2d(Th)((u-0.4*y*(1-y))^2))/sqrt(int2d(Th)((0.4*y*(1-y))^2));
cout <<" L2tau11= " <<L2errortau11<<"L2tau12= "<<L2errortau12<<" L2u="<<L2erroru<< endl;

// Plot each time-step of the 9 components
}
Error:
{
104 : // define DEVSS aux. variable G
105 : solve ComputeG(G11, G12, G21, G22) =
106 : int2d(Th)(G11 * w - dx(u) * w) error operator - <10LinearCombI7MGauche4C_F0E>,
List of choices
( : , )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <14Matrice_CreuseISt7complexIdEE>, <14Matrice_CreuseISt7complexIdEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <14Matrice_CreuseISt7complexIdEE>, <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE>, <14Matrice_CreuseISt7complexIdEE> )
( <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> : <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE>, <NSt7__cxx114listISt5tupleIJSt7complexIdEP13VirtualMatrixIiS3_EbEESaIS7_EEE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <14Matrice_CreuseIdE>, <14Matrice_CreuseIdE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <14Matrice_CreuseIdE>, <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE>, <14Matrice_CreuseIdE> )
( <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> : <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE>, <NSt7__cxx114listISt5tupleIJdP13VirtualMatrixIidEbEESaIS5_EEE> )
( <10FormLinear> : <10FormLinear> )
( <12FormBilinear> : <12FormBilinear> )
( <6C_args> : <6C_args>, <10FormLinear> )
( <6C_args> : <6C_args>, <12FormBilinear> )
( <10LinearCombI7MGauche4C_F0E> : <10LinearCombI7MGauche4C_F0E>, <10LinearCombI7MGauche4C_F0E> )
( <10LinearCombI6MDroit4C_F0E> : <10LinearCombI6MDroit4C_F0E>, <10LinearCombI6MDroit4C_F0E> )
( <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> : <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E>, <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> )
( <10LinearCombI7MGauche4C_F0E> : <10LinearCombI7MGauche4C_F0E> )
( <10LinearCombI6MDroit4C_F0E> : <10LinearCombI6MDroit4C_F0E> )
( <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> : <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E> )
( <4C_F0> : <12TransE_Array>, <12TransE_Array> )
( <4C_F0> : <7E_Array>, <12TransE_Array> )
( <4C_F0> : <7E_Array>, <7E_Array> )
( <4C_F0> : <12TransE_Array>, <7E_Array> )
( <12Add_Mulc_KN_ISt7complexIdEE> : <8Mulc_KN_ISt7complexIdEE>, <3KN_ISt7complexIdEE> )
( <12Add_Mulc_KN_ISt7complexIdEE> : <3KN_ISt7complexIdEE>, <8Mulc_KN_ISt7complexIdEE> )
( <12Add_Mulc_KN_ISt7complexIdEE> : <8Mulc_KN_ISt7complexIdEE>, <8Mulc_KN_ISt7complexIdEE> )
( <8Mulc_KN_ISt7complexIdEE> : <3KN_ISt7complexIdEE> )
( <7Sub_KN_ISt7complexIdEE> : <3KN_ISt7complexIdEE>, <3KN_ISt7complexIdEE> )
( <12Add_Mulc_KN_IdE> : <8Mulc_KN_IdE>, <3KN_IdE> )
( <12Add_Mulc_KN_IdE> : <3KN_IdE>, <8Mulc_KN_IdE> )
( <12Add_Mulc_KN_IdE> : <8Mulc_KN_IdE>, <8Mulc_KN_IdE> )
( <8Mulc_KN_IdE> : <3KN_IdE> )
( <7Sub_KN_IdE> : <3KN_IdE>, <3KN_IdE> )
( <12Add_Mulc_KN_IlE> : <8Mulc_KN_IlE>, <3KN_IlE> )
( <12Add_Mulc_KN_IlE> : <3KN_IlE>, <8Mulc_KN_IlE> )
( <12Add_Mulc_KN_IlE> : <8Mulc_KN_IlE>, <8Mulc_KN_IlE> )
( <8Mulc_KN_IlE> : <3KN_IlE> )
( <7Sub_KN_IlE> : <3KN_IlE>, <3KN_IlE> )
( : )
( : )
( : )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )
( : , )

Error line number 106, in file …edp, before token )

current line = 106
Compile error :
line number :106, )
error Compile error :
line number :106, )
code = 1 mpirank: 0
}

fb77 · February 28, 2025, 5:17pm

I think you can do simply

solve ComputeG11(G11, w) = int2d(Th)(G11 * w) - int2d(Th)(dx(u) * w);
solve ComputeG12(G12, w) = int2d(Th)(G12 * w) - int2d(Th)(dy(u) * w);
solve ComputeG21(G21, w) = int2d(Th)(G21 * w) - int2d(Th)(dx(v) * w);
solve ComputeG22(G22, w) = int2d(Th)(G22 * w) - int2d(Th)(dy(v) * w);

or even simpler by interpolation
G11=dx(u);G12=dy(u);G21=dx(v);G22=dy(v);

Moreover the cross term of the symmetric product should be
(G11 * tau12old + G12 * tau22old + G21 * tau11old + G22 * tau12old))
instead of (G11 * tau12old + G12 * tau22old + G21 * tau12old + G22 * tau11old))

Topic		Replies	Views
Partial periodic boundary condition erroe General Discussion	1	535	June 11, 2020
Compilation error in code " err code 6 , mpirank 0 " application of boundary conditions General Discussion	6	444	January 23, 2023
Code doesn't run on my laptop but do on other's - femlib issues? FreeFEM installation	2	371	October 15, 2021
Assertion fail pointed to DOperator.hpp General Discussion	1	785	January 8, 2020
An error while trying to use periodic boundaries General Discussion	0	338	August 23, 2021

Oldroyd-B model freefem++ code bugs and issues: Assertion fail

Related topics