)
Error line number 103, in file de_kee_stationary_newton.edp, before token )
current line = 103
Compile error :
line number :103, )
error Compile error :
line number :103, )
code = 1 mpirank: 0
```
It looks like there’s a problem multiplying `dvisc * sqrt(...)`. I suspect the problem lies in the fact that the square root is not a linear operation and inside the square root there is my increment `[dux, duy]` to be computed by the solve command. But I would like another opinion on the matter.
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [March 18, 2025, 11:33am UTC](https://community.freefem.org/t/poiseuille-flow-simulation-of-viscoplastic-fluid/3811/2 "2025-03-18T11:33:21Z")
It is difficult to understand such error message, but you are right, the problem is that you must write a variational formulation which is linear with respect to the unknown. For Newton-Raphson, the formula DF(u)(du) = -F(u) is obviously linear with respect to du, hence you must revise your `solve PoiseDKT` to be coherent with the Newton-Raphson algorithm!
Moreover you must separate (i.e. put in different int2d()) terms which are proportional to the unknown and terms which do not depend on the unknown.
---
### Author:  [@eleckbarraki](https://community.freefem.org/u/eleckbarraki)
#### Post date: [April 3, 2025, 3:21pm UTC](https://community.freefem.org/t/poiseuille-flow-simulation-of-viscoplastic-fluid/3811/3 "2025-04-03T15:21:47Z")
Thank you for your answer, I have made the following alterations to the code:
```auto
//
load "iovtk"
//
// Parameters:
real pin = 1.; // input pressure constant
real Re = 0.5; // Reynolds number
real G = 1; // Gamma and theta material parameters of the problem
real th = 15;
real k = th*Re; // DKT number
real Bi = G/(2*Re); // Bingham number
real eps = 1e-4; // term for regularization of viscosity
real err = 0; // error for Newton-Rhapson
real tol = 1e-6;
real n = 20; // number of nodes
int Wall = 1; // pipe wall label
int Inlet = 2; // pipe inlet label
int Outlet = 3; // pipe outlet label
// Boundary
border b1 ( j = 0.0, 1.0 ) { x = 2*j-1; y = -1.0; label=Wall;}
border b2 ( j = 0.0, 1.0 ) { x = 1.0; y = 2*j-1; label=Outlet;}
border b3 ( j = 0.0, 1.0 ) { x = 1-2*j; y = 1.0; label=Wall;}
border b4 ( j = 0.0, 1.0 ) { x = -1.0; y = 1-2*j; label=Inlet;}
// Define the geometry of the problem and the finite element mesh
mesh Mh = buildmesh ( b1(n) + b2(n) + b3(n) + b4(n));
plot(Mh, wait = 1);
// FE space
fespace Vh(Mh,P1); //velocity field
Vh ux,uy,vx,vy,dux,duy;
Vh stress;
Vh gamma,gammae;
Vh visc;
Vh dvisc;
// Macros
macro Grad(ux,uy) [dx(ux), dy(ux) , dx(uy), dy(uy)] // EOM
macro div(ux, uy) (dx(ux) + dy(uy)) // EOM
macro strainrate(ux,uy) [dx(ux),(dy(ux)+dx(uy))/2. , (dy(ux)+dx(uy))/2.,dy(uy)] // EOM
// Set initial conditions and other quantities
ux = 1.-y^4.;
uy = 0.;
dux = 0;
duy = 0;
gamma = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) );
gammae = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) + eps^2.);
visc = Bi/(gammae) + 2.*exp(-2.*k*gamma);
dvisc = Bi/(2.*(gammae)^3.) + 2.*k*exp(-2.*k*gamma)/gammae;
// first an explicit step to feed the Newton loop
solve init([ux,uy],[vx,vy]) =
int2d(Mh)( visc * (strainrate(ux,uy)'*Grad(vx,vy)) - [pin, 0]'*[vx, vy] )
+ on(Wall, ux = 0., uy = 0.) ;
gamma = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) );
gammae = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) + eps^2. );
visc = Bi/gammae + 2.*exp(-2.*k*gamma);
dvisc = Bi/(2.*(gammae)^3.) + 2.*k*exp(-2.*k*gamma)/gammae;
int maxit=50; // maximum number of iterations
for (int i=0;i<=maxit;i++) // Newton-Rhapson DF(u)(du) = -F(u)
{
solve PoiseDKT([dux,duy],[vx,vy]) =
int2d(Mh)(
visc * (strainrate(dux,duy)'*Grad(vx,vy))
- dvisc * (strainrate(ux,uy)'*strainrate(dux,duy)) * (strainrate(ux,uy)'*Grad(vx,vy))
)
- int2d(Mh)( visc * strainrate(ux,uy)'*Grad(vx,vy) - [pin, 0]'*[vx, vy] )
+ on(Wall, dux = 0., duy = 0.) ;
// update u and other quantities
ux[] += dux[]; uy[] += duy[];
gamma = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) );
gammae = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) + eps^2. );
visc = Bi/gammae + 2.*exp(-2.*k*gamma);
dvisc = Bi/(2.*(gammae)^3.) + 2.*k*exp(-2.*k*gamma)/gammae;
err= dux[].linfty + duy[].linfty;
cout << " error = " << err << "\n";
if(err < tol){
cout << " Converges |" <<" iter " << i << " |error = " << err << "\n";
break;
}
if( i=50 && err > 10.){
cout << " Does not converge |" << "error = " << err << "\n";
break;
}
}
plot([u1,u2],cmm="u1,u2 math ",wait=1);
plot (ux, value=1, fill=1, dim=3, wait = 1);
```
I’m trying to run the code but I receive this error message from the terminal:
```auto
86 : int2d(Mh)( visc * (strainrate(ux,uy) [dx(ux),(dy(ux)+dx(uy))/2. , (dy(ux)+dx(uy))/2.,dy(uy)] '*Grad(vx,vy) [dx(vx), dy(vx) , dx(vy), dy(vy)] ) - [pin, 0]'*[vx, vy] ) error operator - <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E>, <10LinearCombI6MDroit4C_F0E>
List of choices
( : , )
.....
```
Then I believe it gives a long list of choices for the elements that can be chosen for the " ’ \* " operator and at the end it confirms that the error is at line 86:
```auto
Error line number 86, in file de_kee_stokes_newton.edp, before token )
current line = 86
Compile error :
line number :86, )
error Compile error :
line number :86, )
code = 1 mpirank: 0
```
To clarify what I’m trying to do, my idea in writing the code is based on the following calculations:
Here D(u) is the symmetric part of the velocity gradient, and u is the velocity. I believe the formulas to be linear w.r.t. the increment du, because both the operator " ’ \* " and the term D(du) are linear. Could the problem be that D(du) contains the derivatives of the increment du? Or am I using the syntax incorrectly somewhere?
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [April 3, 2025, 5:47pm UTC](https://community.freefem.org/t/poiseuille-flow-simulation-of-viscoplastic-fluid/3811/4 "2025-04-03T17:47:11Z")
As I warned you, you must put in different `int2d` terms which are proportional to the unknown and terms which do not depend on the unknown.
Here in `solve init` the term `strainrate(ux,uy)'*Grad(vx,vy)` depends on the unknown, whereas `[pin, 0]'*[vx, vy]` does not. You have to separate them.
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [April 4, 2025, 8:56am UTC](https://community.freefem.org/t/poiseuille-flow-simulation-of-viscoplastic-fluid/3811/5 "2025-04-04T08:56:08Z")
Just a small remark in Newton-Rhapson you can find directly u^{n+1} because the boundary condition on du^{n+1} must be 0 and if you decompose F un 3 part : constant, linear , non nonlinear : F= C + L + NL , when DF(u)(du) = L(du) + DNL(u)(du) and then the Newtow Rhapton become
because we have u^{n+1} = u^{n} +du^{n+1} with DF(u^{n})(du^{n+1}) = -F(u^{n})
DF(u^n)(u^{n+1}) = DF(u^n)(u^{n}) +DF(u^n)(du^{n+1}) = DF(u^n)(u^{n}) -F(u^n)
and with the decomposition you have 🙂
DF(u^n)(u^{n+1}) = DNL(u^n)(u^{n}) -C - NL(u^n)
More easy to code.