# How to set the initial value on the boundary?

**URL:** <https://community.freefem.org/t/how-to-set-the-initial-value-on-the-boundary/2326>\
**Category:** General Discussion\
**Created:** [March 12, 2023, 3:37pm UTC](https://community.freefem.org/t/how-to-set-the-initial-value-on-the-boundary/2326 "2023-03-12T15:37:38Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![cfdyuewu](https://avatars.discourse-cdn.com/v4/letter/c/5f8ce5/32.png) [@cfdyuewu](https://community.freefem.org/u/cfdyuewu)\
**Post date:** [March 12, 2023, 3:37pm UTC](https://community.freefem.org/t/how-to-set-the-initial-value-on-the-boundary/2326/1 "2023-03-12T15:37:38Z")

</div>

I’m working on the example “large fluid problem” by using fem and newton iteration method, I want to set the value of Ux on b6 be 3, because newton method need to solve the incremental value, so I can’t directly assign the dirichlet boundary condition with respect to Ux,Uy,here is the code

load “iovtk”  
verbosity=0;  
// Parameters  
int nn = 15;  
int nnPlus = 5;  
real l = 1.;  
real L = 15.;  
real hSlope = 0.1;  
real H = 6.;  
real h = 0.5;

real beta = 0.01;

real eps = 9.81/303;  
real nu = 1.;  
real numu = nu/sqrt(0.09);  
real nuep = pow(nu,1.5)/4.1;  
real dt = 0.;

real Penalty = 1.e-6;

// Mesh  
border b1(t=0, l){x=t; y=0;}  
border b2(t=0, L-l){x=1.+t; y=-hSlope_t;}  
border b3(t=-hSlope_(L-l), H){x=L; y=t;}  
border b4(t=L, 0){x=t; y=H;}  
border b5(t=H, h){x=0; y=t;}  
border b6(t=h, 0){x=0; y=t;}

mesh Th=buildmesh(b1(nnPlus_nn_l) + b2(nn_sqrt((L-l)^2+(hSlope_(L-l))^2)) + b3(nn\*(H + hSlope\*(L-l))) + b4(nn_L) + b5(nn_(H-h)) + b6(nnPlus_nn_h));  
plot(Th);

// Fespaces  
fespace Vh2(Th, P1b);  
Vh2 Ux, Uy,dux,duy;  
Vh2 Vx, Vy;  
for (int i = 0; i \< Th.nv; ++i)  
{  
if (Th(i).onboundary && Th(i).boundary == 6)  
Ux[i] = 3.0;  
}  
Ux=0+on(b6,Ux=3);  
Vh2 Upx=Ux, Upy=Uy;  
fespace Vh(Th,P1);  
Vh p=0, q,pp;  
Vh Tp, T=35,dT;  
Vh k=log(0.27), kp=k,dk;  
Vh ep=log(0.002104), epp=ep,dep;

fespace V0h(Th,P0);  
V0h muT=1;  
V0h prodk, prode;  
Vh kappa=0.25e-4, stress;  
// Macro  
macro grad(u) [dx(u), dy(u)] //  
macro Grad(U) [grad(U#x), grad(U#y)] //  
macro Div(U) (dx(U#x) + dy(U#y)) //  
macro UgradV(u1,u2,v1,v2) [[u1,u2]‘_[dx(v1),dy(v1)],[u1,u2]'_[dx(v2),dy(v2)]]  
macro ugradv(u1,u2,v) [u1,u2]’\*grad(v);  
// Functions  
func g = (x) \* (1-x) \* 4;  
// Problem  
real alpha = 0.;

// Initialization with stationary solution

plot([Ux, Uy], p, value=true, coef=0.2, cmm=“[Ux, Uy] - p”);  
// Initialization  
dt = 0.05; //probably too big  
int nbiter = 1;  
real coefdt = 0.25^(1./nbiter);  
real coefcut = 0.25^(1./nbiter);  
real cut = 0.01;  
real tol = 0.5;  
real coeftol = 0.5^(1./nbiter);

T = T - 10\*((x\<1)_(y\<0.5) + (x\>=1)_(y+0.1\*(x-1)\<0.5));

// Convergence loop  
real T0 = clock();  
for (int iter = 1; iter \<= nbiter; iter++){  
cout \<\< "Iteration " \<\< iter \<\< " - dt = " \<\< dt \<\< endl;  
alpha = 1/dt;

// Time loop  
real t = 0.;  
for (int i = 0; i \<= 500; i++){  
u\_ite1=Upx;  
u\_ite2=Upy;  
p\_ite=pp;  
T\_ite=Tp;  
k\_ite=kp;  
ep\_ite=epp;

```
  for (int n=0;n<=10;n++){
     solve increment ([du1,du2,dp,dT,dk,dep],[v1,v2,vp,vt,vk,vep])
     =int2d(Th)(alpha*([du1,du2]'*[v1,v2]-dT*vt-dk*vk-dep*vep)
     +0.5*(UgradV(du1,du2,u_ite1,u_ite2))'[v1,v2]
     +0.5*(UgradV(u_ite1,u_ite2,du1,du2))'[v1,v2]
     +0.5*ugradv(v1,v2,dp)
     +(0.09/2)exp(2*k_ite-ep_ite)*(2*k_ite-ep_ite)*(Grad([u_ite1,u_ite2])+Grad([u_ite1,u_ite2])'):Grad([v1,v2])
     +(0.09/2)exp(2*k_ite-ep_ite)*(Grad([du1,du2])+Grad([du1,du2])'):Grad([v1,v2])
     +(9.81/606)*dT*v2-0.5*Div([du1,du2])*vp-0.5*vk*(ugradv(du1,du2,k_ite)+ugradv(u_ite1,u_ite2,dk))-0.5*exp(ep_ite-k_ite)*(dep-dk)*vk+(0.09/2)exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(k_ite),dy(k_ite)]'*[dx(k_ite),dy(k_ite)]*vk
     +(0.09/2)exp(2*k_ite-ep_ite)*2*[dx(k_ite),dy(k_ite)]'*[dx(dk),dy(dk)]*vk
     -(0.09/2)exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(vk),dy(vk)]'*[dx(k_ite),dy(k_ite)]-(0.09/2)exp(2*k_ite-ep_ite)*[dx(vk),dy(vk)]'*[dx(dk),dy(dk)]
     +(0.09/2)*exp(k_ite-ep_ite)*(4*dx(u_ite1)*dx(du1)+4*dy(u_ite2)*dy(du2)+2*(dy(u_ite1)+dx(u_ite2))*(dy(du1)+dx(du2)))*vk
     +(0.09/2)*exp(k_ite-ep_ite)*(dk-dep)*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))*vk
     -0.5*vep*(ugradv(du1,du2,ep_ite)+ugradv(u_ite1,u_ite2,dep))
     -(1.92/2)*exp(-k_ite+ep_ite)*(dep-dk)*vep
     +(0.09/2.6)*exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(ep_ite),dy(ep_ite)]'*[dx(ep_ite),dy(ep_ite)]*vep
     +(0.09/2.6)*exp(2*k_ite-ep_ite)*2*[dx(ep_ite),dy(ep_ite)]'*[dx(dep),dy(dep)]*vep
     -(0.09/2.6)*exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(ep_ite),dy(ep_ite)]'*[dx(vep),dy(vep)]
     -(0.09/2.6)*exp(2*k_ite-ep_ite)*[dx(dep),dy(dep)]'*[dx(vep),dy(vep)]
     +((1.44*0.09)/2)*exp(k_ite-ep_ite)*(4*dx(u_ite1)*dx(du1)+4*dy(u_ite2)*dy(du2)+2*(dy(u_ite1)+dx(u_ite2))*(dy(du1)+dx(du2)))*vep
     +((1.44*0.09)/2)*exp(k_ite-ep_ite)*(dk-dep)*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))*vep
     -0.5*vt*(ugradv(du1,du2,k_ite)+ugradv(u_ite1,u_ite2,dk))
     -0.05*exp(2*k_ite-ep_ite)*(2*dk-dep)*[dx(vt),dy(vt)]'*[dx(T_ite),dy(T_ite)
     -0.05*exp(2*k_ite-ep_ite)*[dx(vt),dy(vt)]'*[dx(dT),dy(dT)
     )
     +int1d(Th,b1,b2,b5)(
        (0.09/2.6)*2.4952*exp(0.5*k_ite)*(0.5*dk)*vep
     )+int1d(Th,b1,b2)(
        0.41*((0.09).^0.25)*0.5*(exp(0.5*k_ite)*(0.5*dk)*u_ite1*v1+exp(0.5*k_ite)*du1*v1)
     )+int2d(Th)(
        alpha*([u_ite1-Upx,u_ite2-Upy]'*[v1,v2]-(T_ite-Tp)*vt-(k_ite-kp)*vk-(ep_ite-epp)*vep)
        +0.5*(UgradV(u_ite1,u_ite2,u_ite1,u_ite2))'*[v1,v2]+0.5*(UgradV(Upx,Upy,Upx,Upy))'*[v1,v2]
        +0.5*grad(p_ite+pp)'[v1,v2]
        +0.09*0.5*exp(2*k_ite-ep_ite)*(Grad([u_ite1,u_ite2])+(Grad([u_ite1,u_ite2]))'):Grad([v1,v2])
        +0.09*0.5*exp(2*kp-epp)*(Grad([Upy,Upy])+(Grad([Upx,Upy]))'):Grad([v1,v2])
        +(9.81/606)*(T_ite+Tp-70)*v2-0.5*Div([u_ite1+Upx,u_ite2+Upy])*vp
        -0.5*(ugradv(u_ite1,u_ite2,k_ite))*vk-0.5*(ugradv(Upx,Upy,kp))*vk
        -0.5*exp(ep_ite-k_ite)*vk-0.5*exp(epp-kp)*vk
        +(0.09/2)*vk*(exp(2*k_ite-ep_ite)*grad(k_ite)'*grad(k_ite))
        +(0.09/2)*vk*(exp(2*kp-epp)*grad(kp)'*grad(kp))
        -(0.09/2)*(exp(2*k_ite-ep_ite)*grad(k_ite)'*grad(vk))
        -(0.09/2)*(exp(2*kp-epp)*grad(kp)'*grad(vk))
        +(0.09/2)*exp(k_ite-ep_ite)*vk*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))
        +(0.09/2)*exp(kp-epp)*vk*(2*dx(Upx)*dx(Upx)+2*dy(Upy)*dy(Upy)+(dy(Upx)+dx(Upy))*(dy(Upx)+dx(Upy)))
        -0.5*(ugradv(u_ite1,u_ite2,ep_ite))*vep-0.5*(ugradv(Upx,Upy,epp))*vep
        -1.92*0.5*exp(-k_ite+ep_ite)*vep--1.92*0.5*exp(-kp+epp)*vep
        +(0.09/2.6)*vep*(exp(2*k_ite-ep_ite)*grad(ep_ite)'*grad(ep_ite))+(0.09/2.6)*vep*(exp(2*kp-epp)*grad(epp)'*grad(epp))
        -(0.09/2.6)*(exp(2*k_ite-ep_ite)*grad(ep_ite)'*grad(vep))-(0.09/2.6)*(exp(2*kp-epp)*grad(epp)'*grad(vep))
        +(1.44*0.09/2)*exp(k_ite-ep_ite)*vep*(2*dx(u_ite1)*dx(u_ite1)+2*dy(u_ite2)*dy(u_ite2)+(dy(u_ite1)+dx(u_ite2))*(dy(u_ite1)+dx(u_ite2)))
        +(1.44*0.09/2)*exp(kp-epp)*vep*(2*dx(Upx)*dx(Upx)+2*dy(Upy)*dy(Upy)+(dy(Upx)+dx(Upy))*(dy(Upx)+dx(Upy)))
        -0.5*(ugradv(u_ite1,u_ite2,T_ite))*vt-0.5*(ugradv(Upx,Upy,Tp))*vt
        -0.05*(exp(2*k_ite-ep_ite)*grad(T_ite)'*grad(vt))-0.05*(exp(2*kp-epp)*grad(Tp)'*grad(vt))
     )+int1d(Th,b1,b2,b5)(
        (0.09/2.6)*2.4952*exp(0.5*k_ite)*vep+(0.09/2.6)*2.4952*exp(0.5*kp)*vep)
     +int1d(Th,b1,b2)(0.41*((0.09).^0.25)*exp(0.5*k_ite)*u_ite1*v1*0.5+0.41*((0.09).^0.25)*exp(0.5*kp)*Upx*v1*0.5)
     +on(b5,du1=0,du2=0)
     +on(b6,du1=0,du2=0,dep=0,dk=0,dT=0)+on(b2, du2=-du1*N.x/N.y,dT=0)+on(b1,du2=0,dT=0);
     u_ite1+=du1;
     u_ite2+=du2;
     p_ite+=dp;
     T_ite+=dT;
     k_ite+=dk;
     ep_ite+=dep;
  }
  Upx=u_ite1;
  Upy=u_ite2;
  pp=p_ite;
  Tp=T_ite;
  kp=k_ite;
  epp=ep_ite;
  plot([Upx, Upy], p, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);

```

}

// Check  
if (iter \>= nbiter) break;  
}  
cout \<\< "Total Time = " \<\< clock()-T0 \<\< endl;  
how to do this?