# Non homogeneous boundary and initial conditions

**URL:** https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154
**Category:** General Discussion
**Created:** [October 12, 2019, 9:10am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154 "2019-10-12T09:10:13Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![Aicha](https://avatars.discourse-cdn.com/v4/letter/a/ec9cab/32.png) [@Aicha](https://community.freefem.org/u/Aicha)
#### Post date: [October 12, 2019, 9:10am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/1 "2019-10-12T09:10:13Z")

</div>

Hi,  
i try to write a code who calculate the solution of the problem:

\partial\_t u + u \partial\_x u - \nu \partial^2\_x u =1+x, \ x \in [0,1], t \in [0,1]

with initial condition:

u(x,0)= -4 \nu \dfrac{x+3/2}{ (x+3/2)^2+\nu } - (x+2)

and boundary conditions:

u(0,t)= -4 \nu \dfrac{1+1/2 e^t}{ (1+1/2 e^t)^2 +\nu }-(1+e^t),

u(1,t)= -4 \nu \dfrac{ 2+1/2 e^t }{ (2+1/2 e^t)^2 +\nu} -(2+e^t).

The exact solution of this problem is:

u(x,t)=-4 \nu \dfrac{1+x+1/2 e^t} { (1+x+1/2 e^t)^2+\nu } -(1+x+e^t)

We can take \nu=1..

The goal is to calculate the approximate solution with FF+ and compare it to the exact solution compring the graphs of the two functions. I try to write the Following code but there is an error and i Don’t find how correct it.  
Thank you in advance to the help.

```auto
load "iovtk"

verbosity=1;
real nu=1.;
int n = 100; 
real a=0, b=1;
mesh Th=square(n,1,[a+x*(b-a),y/10]);
//plot(Th,wait=false);
fespace Vh(Th,P2,periodic=[[1,x],[3,x]]);
Vh u,v;
u = 0;
real T=0;
func f = 1+x;

func real uex(real t) //The exact solution
{
return -4*nu* (1.+x+(1./2.)*exp(t))/ ( (1.+x+(1./2.)*exp(t))^2+nu) - (1.+x+exp(t));
};

func up= -4*nu* (x+(3./2.))/((x+(3./2.))^2+nu)-(x+2.);//Initial condition

func real g(real t) //Boundary condition u(x=0,t)
{
return -4*nu*(1.+(1./2.)*exp(t))/((1+(1./2.)*exp(-t))^2+nu) -(1+exp(t));
};

func real r(real t) //Boundary condition u(x=1,t)
{
return -4*nu*(2.+(1./2.)*exp(t))/((2.+(1./2.)*exp(t))^2+nu) -(2+exp(t));
};

Vh fh=f;
Vh gh=g;
Vh rh=r;
Vh uph=up; 
real dt = 0.1;
Vh g0,g1;
Vh nith = 0.0; // for nonlinear iterations

problem Pb(u,v) = int2d(Th)( u*v/dt + nu*dx(u)*dx(v) ) - int2d(Th)( up*v/dt + f*v)
+ int2d(Th) (g0*v) + int2d(Th) (g1*u*v) + int2d(Th) (nith*dx(u)*v)+ on(2,u=g)+ on(4,u=r);

int k=0;        
int noit = 10; // maximum of iterations
real errnlpb = 0; // error of iteration of Pb
real normuh = 0;
 
ofstream fout("uh.txt");
fout << "# t uh" << endl;

ofstream fouti("U.txt");
fouti << "# t U" << endl;

for(int it=0; it<1.; ++it)
{
     T += dt; 
    // Nonlinear iterations for Pb solution
    nith = u; // start nonlinear iterations with last solution for Pb
    g0=u*dx(u);
    g1=dx(u);
    cout << " Solving nonlinear problem.\n";
    for(k=1;k<noit; ++k) {
     up[]=u[];
     Pb;
      normuh = sqrt(int1d(Th,1)(square(u))); // L^2 nor of the solution of Pb
      errnlpb = sqrt(int1d(Th,1)(square(u - nith))); // L^2 norm of the error in nonlinear iteration of Pb
      cout << " " << k << ": error = " << errnlpb << endl;
      if(errnlpb <= 1E-5*normuh) // stopping criteria. If relative error
      break; // is less than 1E-5 stop nonlinear iterations of Pb
      nith = u; // prepare new iteration
      g0=u*dx(u);
      g1=dx(u);
      }
    // Check convergence for Pb  
    if(k >= noit) {
      cout << "Nonlinear solver for Pb did not converged.\n";
      cout << "Error = " << errnlpb << endl;
      cout << "|| u|| = " << normuh << endl;
      exit(1); // exit the program with error code 1
      }
Vh uh=uex;

 plot(u,wait=1,dim=1,fill=true,cmm=it,value=1);

 plot(uh,wait=1,dim=1,fill=true,cmm=it,value=1);

fout.flush;

fouti.flush;

}

```

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [October 20, 2019, 3:35pm UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/2 "2019-10-20T15:35:05Z")

</div>

You are solving a time-dependant nonlinear problem. There are tools out there to help you so that you don’t have to reimplement everything yourself. You should use TSSolve, cf. this linear [example](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/heat-2d-TS-RHS-PETSc.edp) + this [paper](https://arxiv.org/pdf/1806.01437.pdf). You need to supply two FreeFEM functions to update the Jacobian of your problem as well as evaluate the residual.

---

<div class="post-metadata">

### Author: ![Aicha](https://avatars.discourse-cdn.com/v4/letter/a/ec9cab/32.png) [@Aicha](https://community.freefem.org/u/Aicha)
#### Post date: [November 1, 2019, 10:19am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/3 "2019-11-01T10:19:48Z")

</div>

Hi  
i try the idea but it’s the same, i obtain no good results. Can any one help me please to write a code who resolve the problem  
$$\partial\_t u + u \partial\_x u - \nu \partial\_x^2 u = 1+x, \ t \>0, x \in ]0,1[,$$  
with initial condition  
u(x,0)= \dfrac{- 2 \nu}{x+3/2} +2+x,  
and boundary conditions  
u(0,t)= \dfrac{-2 \nu}{1+1/2 e^{-t}} +1+e^{-t} and u(1,t)= \dfrac{-2 \nu}{2+1/2 e^{-t}} +2+e^{-t}.

and compare the approximate solution with the exact solution  
$$u(x,t)= -2 \nu \dfrac{1}{1+x+1/2 e^{-t}} +x+1+e^{-t}$$

Thank’s in advance to the help.

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 11, 2023, 8:29am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/4 "2023-05-11T08:29:27Z")

</div>

Here is an example code I created for TSSolve with Navier Stokes that may be helpful for you (or anyone else). It is not at all optimal, but does show the basic steps involved using TSSolve for nonlinear systems.  
[NS-TSSolve.edp](https://community.freefem.org/uploads/short-url/xCH846sE8CUDGD9r8MdNpY7spOW.edp) (2.9 KB)

---

<div class="post-metadata">

### Author: ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [May 11, 2023, 10:49am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/5 "2023-05-11T10:49:48Z")

</div>

How does this compare to SNESSolve? The one example of this I thought  
was useful and it has been producing believable ( I have not checked details )  
results for me in a simple drift-diffusion problem.

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 11, 2023, 11:00am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/6 "2023-05-11T11:00:57Z")

</div>

TSSolve uses SNESSolve to resolve the nonlinear problems at each time step.

---

<div class="post-metadata">

### Author: ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [May 11, 2023, 11:04am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/7 "2023-05-11T11:04:48Z")

</div>

Since the range of t appears to be limited and the mesh is not likely to change  
unpredictably, with time, what is wrong with just using the y axis for time?  
Since the OP provides an analytical solution, I’ll see if I can set it up this  
way 🙂 SNESSolve shouldn’t know the difference between y and t…

---

<div class="post-metadata">

### Author: ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [May 11, 2023, 1:06pm UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/8 "2023-05-11T13:06:20Z")

</div>

Frankly this looks too good to be accurate but the reason I’m responding is  
to test my knowledge and scripts 🙂 See if the attached is helpful. The relevant function  
and jacobian are included but a lot of the snes stuff is in my macro which I can post  
but would like to clean up first. Fortuitous or accidental agreements may be quite  
common and there are likely errors here so comments appreciated. Thanks.  
Note that “y” is time and the diff is the analytical soln minus the FF soln.  
See the code for details.

```auto
l2=4.354910e-09 l2soln=2.196520e+01 mssoln=2.178491e+01 sd-1.802880e-01
times: compile 1.298770e-01s, execution 1.074365e+00s, mpirank:0
 ######## We forget of deleting 27 Nb pointer, 0Bytes , mpirank 0, memory leak =508272
 CodeAlloc : nb ptr 8250, size :716256 mpirank: 0
Ok: Normal End

```

 ![forum11soln](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/9/96d1506875a5255d43b9f5f1af1a72d78e3b47ea.jpeg)  
 ![forum11](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/dc0c0b576ca148ad4aaaaf24e3ff8fe4b3c0ade2.jpeg)  
[forum11snes.edp](https://community.freefem.org/uploads/short-url/kvDn6xxLYwn7FHBn2SjoA5LlEjH.edp) (2.4 KB)

```auto

// forum11snes.edp
// https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/7

load "PETSc"
load "MUMPS_mpi"
 macro dimension()2// EOM
include "macro_ddm.idp"
include "cube.idp"
include "mysness.idp"

// needed for mysness.idp
include "myboxes.idp"
include "surfaces.idp"
include "mjmio.idp"

func Pk = P1; // Finite-element space

int nx = getARGV("-nx" , 20) ; //
int nt = getARGV("-nt" , 20) ; //
real nu = getARGV("-nu" ,1.0) ; //
int split= getARGV("-split", 1);

// b,r,t,l
int[int] lbl=[1,2,3,4];
mesh Th = square(nx,nt,label=lbl);

// boundary funcs and known solutions
func real exact(real xx, real t) { 
real et=exp(t);
real p1=1+xx+.5*et;
return -4.0*nu*(p1)/( p1*p1+nu)-(1+xx+et);
} // exact
func real uinit(real xx)
	{real p1=xx+1.5; return -4.0*nu*(p1)/(p1*p1+nu)-(xx+2); }
func real uzt(real t){real et=exp(t); real p1=1+.5*et; 
	return -4.0*nu*(p1)/(p1*p1+nu)-(1+et); }
func real uonet(real t){real et=exp(t); real p1=2+.5*et; 
	return -4.0*nu*(p1)/(p1*p1+nu)-(2+et); }

fespace Vh(Th,Pk);
// a bunch of stuff for my macros I need to work out yet. 
// this defines u and A etc
Mat A;
buildMat(Th, split, A, Pk, mpiCommWorld);
int solnguessidx=0;

Vh u = 0;
Vh soln=exact(x,y);
mesh ThBackup=Th;

// the boundary conditions 
macro EQUBC +on(1,phi=uinit(x)) +on(2,phi=uonet(y)) +on(4,phi=uzt(y)) // EOM

// needed functions, a bit idiosyncratic maybe need to fix macro
// but should be more or less decipherable
varf equpro(phi,vphi)
=int2d(Th)(dy(phi)*vphi+u*dx(phi)*vphi+nu*dx(phi)*dx(vphi))
-int2d(Th)(vphi*(1+x))
EQUBC
;
// use u here not phi... 
varf equres(phi,vphi)
=int2d(Th)(dy(u)*vphi+u*dx(u)*vphi+nu*dx(u)*dx(vphi))
-int2d(Th)(vphi*(1+x))
EQUBC ;
// diff wrt "u"
varf equj(phi,vphi)
=int2d(Th)(dy(phi)*vphi+dx(phi)*vphi+nu*dx(phi)*dx(vphi))
EQUBC ;

snesjac(equJ,equj)
// define the evaluate function
sneseval(equRes,equres)
// the third parameter is probably not right in general
mysnesflags(u,equpro,equJ,equRes,equj,0);
// my macro returns u;
Vh diff=soln-u;
real l2=int2d(Th)(diff*diff);
real l2soln=int2d(Th)(soln*soln);
real mssoln=int2d(Th)(soln);
mssoln=mssoln*mssoln;
cout<<" l2="<<l2<<" l2soln="<<l2soln<<" mssoln="<<mssoln<<" sd-"<<(l2soln-mssoln)<<endl; cout.flush;

plot(soln,cmm="analyical", fill=1,wait=1,value=1);
plot(u,cmm="ff soln",fill=1,wait=1,value=1);
plot(diff,cmm="diff", fill=1,wait=1,value=1);
mmsave("forummacrotest.txt",u,Vh);

```

---

<div class="post-metadata">

### Author: ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [May 11, 2023, 1:09pm UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/9 "2023-05-11T13:09:35Z")

</div>

here is the output file fwiw,

[forummacrotest.zip](https://community.freefem.org/uploads/short-url/pByOWBN5FYe8jTQd8lsDo05EdzZ.zip) (2.4 KB)

---

<div class="post-metadata">

### Author: ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)
#### Post date: [May 6, 2024, 8:46am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/10 "2024-05-06T08:46:39Z")

</div>

Hi cmd,  
thanks for the example. I tried to understand this code and I ran the code to see the result, but, it looks weird and can’t converge.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/a/a12fed06e6f67b48503823b6d2fdc572b9e6208f.png)

I don’t know what’s wrong, do you have any comments?

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 6, 2024, 10:31am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/11 "2024-05-06T10:31:56Z")

</div>

That is just an unhelpful warning that appears in parallel FF codes. It is not relevant here. You should run it with `-v 0`.

Mind the partition of unity!

---

<div class="post-metadata">

### Author: ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)
#### Post date: [May 6, 2024, 10:42am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/12 "2024-05-06T10:42:48Z")

</div>

Thanks for your quick reply, I use -v 0, and it works like this:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/549439222be319958039687dbacdc1df58b2e3d9.png)

it stuck on step 2.  
and if run with -n 1, there no such problem.

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 6, 2024, 11:04am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/13 "2024-05-06T11:04:46Z")

</div>

Sorry. This was my mistake. Here’s a code that works.  
[NS-TSSolve(1).edp](https://community.freefem.org/uploads/short-url/7d0ZQv9ceaGIKoyvXEig3lRLBwz.edp) (3.0 KB)

---

<div class="post-metadata">

### Author: ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)
#### Post date: [May 6, 2024, 11:56am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/14 "2024-05-06T11:56:16Z")

</div>

Thanks for the fix, it works well now.  
What do the parameters `inverse` and `exchange` mean?  
In this example, the `funcF` and `funcJ` uses ` ChangeNumbering(TS, ub[], in, inverse = true, exchange = true);`  
I know that `ub[]` is a FreeFEM array, and `in` is a PETSc numbering, what happened to them after the changenumbering?

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 6, 2024, 12:16pm UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/15 "2024-05-06T12:16:24Z")

</div>

`ChangeNumbering` handles the correspondence between the FreeFEM and PETSc numbering schemes.  
The `exchange` argument updates ghost values if true.  
The `inverse` argument controls whether the renumbering is from FreeFEM to PETSc (false) or PETSc to FreeFEM (true).

You can find more details [here](https://joliv.et/FreeFem-tutorial/main.html).

---

<div class="post-metadata">

### Author: ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)
#### Post date: [May 6, 2024, 12:19pm UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/16 "2024-05-06T12:19:24Z")

</div>

Thanks, it is more clear now.  
Another problem, in each step, if I want to solve fluid system(on Fluid mesh) first, and then solve the thermal field(Fluid+Solid), could it be implemented in TSSolve?

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 6, 2024, 12:28pm UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/17 "2024-05-06T12:28:21Z")

</div>

Yes. The easiest way (implementation-wise) will be to solve that whole coupled system in one shot. You should be able to get better performance if you do a splitting like you described. In that case, you could just time-step the NSE using a code very similar to the example, but add an additional routine to compute the thermal field at each time step in the function for the residual provided to `TSSolve`.

---

<div class="post-metadata">

### Author: ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)
#### Post date: [May 7, 2024, 10:30am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/18 "2024-05-07T10:30:49Z")

</div>

Hi cmd,  
I tried to change the boundary condition to  
`+ on(3, dub =ub- 1, duby = uby) + on(1,2,4, dub = ub, duby = uby);`  
and  
`+ on(3, ub = 1, uby = 0) + on(1,2,4, ub = 0, uby = 0);`for Stokes  
to compute a lid-driven cavity flow. I think it would work as previous, however, the TSSolve can’t converge for this modify.I tried other -ts\_type, but failed.  
After that, I added `+epsp*dubp*vp` to vJ where epsp is 1e-6. It works but I am not sure. Is this a right way?

---

<div class="post-metadata">

### Author: ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [May 7, 2024, 10:48am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/19 "2024-05-07T10:48:36Z")

</div>

Yes this is expected because the pressure is only defined up to a constant without an outflow condition or other constraint on the pressure.

The penalization you describe is probably the easiest solution, though it is not the best from the perspective of actually conserving mass. I think a better approach is to introduce a scalar Lagrange multiplier to constrain the average pressure to a fixed value. This requires you to augment the Jacobian matrix with a single row and column describing the constraint.

Are you actually interested in the cavity flow or was this just a test case?

---

<div class="post-metadata">

### Author: ![quentin](https://avatars.discourse-cdn.com/v4/letter/q/e47c2d/32.png) [@quentin](https://community.freefem.org/u/quentin)
#### Post date: [May 7, 2024, 11:14am UTC](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154/20 "2024-05-07T11:14:12Z")

</div>

Just a test case, indeed, I am working on a convection problem with a solid phase participating in the heat transfer, and first step is to run DNS with multiple parameters to get data for analysis. I have done my serial code based on Hecht’s code, but too slow on high Rayleigh number, so I am learning to parallelize it.

[Next page](https://community.freefem.org/t/non-homogeneous-boundary-and-initial-conditions/154.md?page=2)