# Getting Problem to compute Intermediate velocity of Chan-Hilliard-Navier-stokes equation

**URL:** <https://community.freefem.org/t/getting-problem-to-compute-intermediate-velocity-of-chan-hilliard-navier-stokes-equation/3646>\
**Category:** General Discussion\
**Created:** [December 18, 2024, 5:52am UTC](https://community.freefem.org/t/getting-problem-to-compute-intermediate-velocity-of-chan-hilliard-navier-stokes-equation/3646 "2024-12-18T05:52:02Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Monirul25](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/monirul25/32/3287_2.png) [@Monirul25](https://community.freefem.org/u/Monirul25)\
**Post date:** [December 18, 2024, 5:52am UTC](https://community.freefem.org/t/getting-problem-to-compute-intermediate-velocity-of-chan-hilliard-navier-stokes-equation/3646/1 "2024-12-18T05:52:02Z")

</div>

Dear All, i am trying to solve Chan-Hilliard-Navier-Stokes equation by Discontinuous Galerkin Pressure Correction method. For that i have to compute an intermediate velocity vtilde. I am getting problem to solve the equation.

Here, is my equation

 ![eqn](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/d14f1aa7f2cea30c7aae760053c64b1cedd5fb51.png)  
Here, are the bilinear forms:  
 ![b_one](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c26381dd505488a83aaa5bf3cdf067e22a3c34a6.png)  
 ![b_two](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/a/ae88ddf32b975603358888091d4150ffba8a4bbc.png)  
 ![b_three](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/7/7dad1184c4da323521ea66f1c5f4b4e6576f1105.png)  
 ![b_four](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/1/152fa60bc58833d3a4d7562d18a3ef27fc59d589.png)  
Here, is my FF++ code

> Blockquote  
> Given that v1=0,v2=0 on boundary;  
> so that v1tilde=0, v2tilde=0 on boundary;

real nu=1;// fluid viscosity  
real b=1./2;

// To solve intermediate velocity  
varf pb4tilde(v1tilde,v2tilde,vv1,vv2)=int2d(Th)(v1tilde_vv1/dt+v2tilde_vv2/dt)

```
	                                 // Convection terms
                                        +int2d(Th)(v1tilde0*dx(v1tilde)*vv1+v2tilde0*dy(v1tilde)*vv1+v1tilde0*dx(v2tilde)*vv2+v2tilde0*dy(v2tilde)*vv2) 
                                        +intalledges(Th)(((nTonEdge-1)*real(mean(v1tilde0)*N.x+mean(v2tilde0)*N.y<0.)*abs(N.x*mean(v1tilde0)+N.y*mean(v2tilde0))*(-jump(v1tilde)*vv1-jump(v2tilde)*vv2))/nTonEdge)
                                        +int1d(Th)(real(v1tilde0*N.x+v2tilde0*N.y<0.)*abs(N.x*v1tilde0+N.y*v2tilde0)*(v1tilde*vv1+v2tilde*vv2)) // inflow boundary
										
										+int2d(Th)(((dx(v1tilde0)+dy(v2tilde0))*v1tilde*vv1+(dx(v1tilde0)+dy(v2tilde0))*v2tilde*vv2)/2.)
										-intalledges(Th)(b*(jump(N.x*v1tilde0+N.y*v2tilde0)*(mean(v1tilde)*mean(vv1)+(jump(v1tilde)*jump(vv1))/4.)   
                                        +jump(N.x*v1tilde0+N.y*v2tilde0)*(mean(v2tilde)*mean(vv2)+(jump(v2tilde)*jump(vv2))/4.))/nTonEdge*(nTonEdge-1)) 
                                        -int1d(Th)(b*((N.x*v1tilde0+N.y*v2tilde0)*(v1tilde*vv1+v2tilde*vv2)))
										
										
                                        +int2d(Th)(nu*(dx(v1tilde)*dx(vv1)+ dy(v1tilde)*dy(vv1)+dx(v2tilde)*dx(vv2)+ dy(v2tilde)*dy(vv2)))
									    -int2d(Th)(p0*(dx(vv1)+dy(vv2)))+int2d(Th)(u0*(dx(w)*vv1+dy(w)*vv2))
										+intalledges(Th)(
                                            // loop on all edges of all triangles
                                             // the edges are seen nTonEdge times so we divide by nTonEdge
                                            // remark: nTonEdge = 1 on border edges and = 2 on internal edges
                                           // in a triangle, the normal is the exterior normal
                                          // def: jump = external - internal value; on border, external value = 0
                                          // mean = (external + internal value)/2, on border just internal value

                                          // Main computation for edge integrals
                                       ((nu*(-mean(dn(v1tilde))*jump(vv1)-mean(dn(vv1))*jump(v1tilde)+(1000./lenEdge)*jump(v1tilde)*jump(vv1))
                                        +nu*(-mean(dn(v2tilde))*jump(vv2)-mean(dn(vv2))*jump(v2tilde)+(1000./lenEdge)*jump(v2tilde)*jump(vv2))))/nTonEdge*(nTonEdge - 1))
										-int1d(Th)(nu*dn(v1tilde)*vv1)-int1d(Th)(nu*dn(v2tilde)*vv2) //-int1d(Th)(nu*dn(vv1)*v1tilde)-int1d(Th)(nu*dn(vv2)*v2tilde)
								        +intalledges(Th)((mean(p0)*jump(N.x*vv1+ N.y*vv2)-mean(u0)*jump(w)*mean(N.x*vv1+ N.y*vv2))/nTonEdge*(nTonEdge-1))
										+int1d(Th)(p0*(N.x*vv1+N.y*vv2));
										
										
										
 solve utilde(v1tilde,v2tilde,vv1,vv2,solver=UMFPACK) = pb4tilde
                                                      // RHS terms
							                          -int2d(Th)(v1tilde0*vv1/dt+g1(t+dt/2.)*vv1)
                                                      -int2d(Th)(v2tilde0*vv2/dt+g2(t+dt/2.)*vv2);
													  
													  
// Solve utilde
utilde;	// Optimzed Error comming here

```

> Blockquote

Here, is the error

 ![Error](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/0ba5daa4f7f68fc7411aca2a86102bf0533aa7cd.png)

**As per my understanding, it is coming as the matrix i not invertible but i am done exactly what same with exiting equations.**

Can you help me please to solve this error.

\*\*Here, \mu^n\_h= w^n\_h , is computed already computed from previous steps. Also, c^{n-1}\_h= c^{o} is known from previous step. \*\*  
**Code is running and giving good results if i change the space of velocity Xh(Th, P2) (Continuous space) but optimized error coming by considering vtilde in Vh(Th, P2dc).**