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

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


Here, are the bilinear forms:




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)(v1tildevv1/dt+v2tildevv2/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

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).