Updating a matrix for a transient convection problem without using convect()

General Discussion
1

I am having difficulty updating the LHS of my problem of the form

L(u,v) + N(\bar{u},u,v) = f

where L is linear, N is not such that \bar{u} = F(u) (for some F) and is linearized by either iterating or otherwise. (In this case below, it’s just the navier stokes equations in the usual skew-symmetric formulation of the advection. Time-stepping with backwards-euler, I’m simply linearizing N with the previous time-step (as an example.) I’ve tried several things, but I can’t seem to find a clear way to update the system to reflect a changing \bar{u}- I feel like I’m probably making a silly mistake somewhere.

EDIT: This code should work fine now, as written. New code below highlights the same problem I stated above, essentially, but it has to do with the Mat-type wrappers used in the parallel sparse solvers.

load "iovtk"

int nn=50;

mesh Th=square(nn,nn);

macro GradGrad(u,v) (dx(u)*dx(v) + dy(u)*dy(v)) //
macro Convection(w1,w2,u,v) (v*(w1*dx(u) + w2*dy(u)) - u*(w1*dx(v) + w2*dy(v)))//

fespace Xh(Th,P2);
fespace Qh(Th,P1);

int DOFs = 2*Xh.ndof + Qh.ndof;
int vDOF = Xh.ndof;
int pDOF = Qh.ndof;

Xh u1, u2, v1, v2;
Xh u1m, u2m;
Qh p,q;

func lid = (4.0*x*(1-x))^0.01;

real TMax = 5;
int Nsteps = 5000;
real dt = TMax/Nsteps;

real nu = 1.0/4000;

//Initialize the data
u1 = 0;
u2 = 0;
u1m = 0;
u2m = 0;

Xh u,v;
varf varfNSEU1V1(u,v) = int2d(Th)(u*v + nu*dt*GradGrad(u,v) + dt*Convection(u1m,u2m,u,v)) + on(1, u=lid) + on(2,3,4, u=0);
varf varfNSEU2V2(u,v) = int2d(Th)(u*v + nu*dt*GradGrad(u,v) + dt*Convection(u1m,u2m,u,v)) + on(1,2,3,4, u=0);
varf varfNull(u,v)   = int2d(Th)(u*v + nu*dt*GradGrad(u,v) + dt*Convection(u1m,u2m,u,v));
varf varfNSERHS(u,v) = int2d(Th)(u*v);



varf varfqdiv1(u,q) = -int2d(Th)(dt*q*dx(u)); 
varf varfqdiv2(u,q) = -int2d(Th)(dt*q*dy(u));
varf varfPmass(p,q)  = int2d(Th)(1e-10*dt*p*q);

matrix pMass = varfPmass(Qh,Qh);
matrix Bx    = varfqdiv1(Xh,Qh);
matrix By    = varfqdiv2(Xh,Qh);
matrix velmass = varfNSERHS(Xh,Xh);
matrix RHSPROJ = [[velmass, 0, 0],[0, velmass,0],[0,0,pMass]];



	 

real[int] b0(DOFs); //NSE * xx = b
real[int] b(DOFs);



real[int] xx(DOFs);


//RHS Terms
real[int] bu1 = varfNSEU1V1(0,Xh);
real[int] bu2 = varfNSEU2V2(0,Xh);
//Labels (1 or 0)
real[int] bu1bdy(vDOF);
real[int] bu2bdy(vDOF);
real[int] bu3bdy(vDOF);

for(int ii=0; ii<vDOF; ii++){
   b0[ii]      = bu1[ii];
   b0[ii+vDOF] = bu2[ii];
}
b = 1.0*b0;

for(int ii=0; ii<pDOF; ii++){
	b[ii + 2*vDOF] =0;
}

for(int ii=0; ii<vDOF; ii++){
	if(bu1[ii] == 0){
		bu1bdy[ii] = 1;
	}
	if(abs(bu1[ii])>0){
		bu1bdy[ii] = 0;
	};
	if(bu2[ii] == 0){
		bu2bdy[ii] = 1;
	}
	if(abs(bu2[ii])>0){
		bu2bdy[ii] = 0;
	};
}

//Solve

matrix U1V1, U2V2, UVN, NSE, Null;
int k = 0;



for(int nt = 0; nt<Nsteps; nt++){
//This is where I'm attempting to update the matrix NSE
	 varf varfNSEU1V1(u,v) = int2d(Th)(u*v + nu*dt*GradGrad(u,v) + dt*Convection(u1m,u2m,u,v)) + on(1, u=lid) + on(2,3,4, u=0);
	 varf varfNSEU2V2(u,v) = int2d(Th)(u*v + nu*dt*GradGrad(u,v) + dt*Convection(u1m,u2m,u,v)) + on(1,2,3,4, u=0);
	 varf varfNull(u,v) = int2d(Th)(u*v + nu*dt*GradGrad(u,v) + dt*Convection(u1m,u2m,u,v));
	 
	 U1V1  = varfNSEU1V1(Xh,Xh);
	 U2V2  = varfNSEU2V2(Xh,Xh);
	 UVN   = varfNull(Xh,Xh);
	
	 NSE = [[U1V1, 0, Bx'],
			[0, U2V2, By'],
			[Bx, By, pMass]];	

	 Null = [[UVN, 0, Bx'],
			[0, UVN, By'],
			[Bx, By, pMass]];	
			
	set(NSE,solver=sparsesolver);
	xx = NSE^-1*b;
	for(int ii=0; ii<vDOF; ii++){
		u1m[][ii] = xx[ii];
		u2m[][ii] = xx[vDOF+ii];
		}
	for(int ii=0; ii<pDOF; ii++){
		p[][ii] = xx[2*vDOF + ii];
	}
	//Corrected RHS Update
	b = RHSPROJ*xx;
	b += b0;
	
	if(mpirank == 0){
	 for(int ii = 0; ii < 10; ii++){
	 k = nt%50;
	 if(k==0){
		savevtk("nse_test_" + nt + ".vtu",Th,[u1m,u2m,0],p,dataname="Velocity Pressure");
	 }
	 cout << xx[ii]<<endl;
	 cout << DOFs << endl;
	 }
	 }
}

UPDATE- I managed to fix the issue, as stated above. The matrix objects can and are updated just fine as I’ve done, I just made a mistake when re-assigning my RHS. I’ve corrected the above code to function as I intended.

However, I wrote this question in lieu of my full question, assuming I had honed in on my actual problem. When the problems become too large for UMFPACK, I can use GMRES by including the following

Mat matNSE(NSE);
set(matNSE,'settings..')
xx = matNSE^-1 * b;

If I loop, and update the NSE matrix, I can’t seem to update the matNSE without getting an error.

Help? (haha)