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

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

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)