Exporting mass matrix and load vector?

[Jared Barber]

Hello,

I was trying to export mass matrices and load vectors when I saw something I did not expect.

Looking at the dbg file it appears the mass matrix for a given mesh refinement corresponds to the heading "SPACE JACOBIAN".

After the mass matrix appears a set of RHS and SOLUTION vectors. It seems that for a linear problem (like the one below that I used to check this out) there are nested iterations with the solution of the outermost iteration printed out at each outermost iteration step. Hence we see a handful of rhs's and solution's for each grid.

The issue is that the rhs, which I thought was the load vector, seems to change during these outermost iterations. I thought that for a given grid the load vector for a linear problem (and the below problem is linear) should not change.

Hopefully someone can help me understand what's going on here. Is RHS the load vector? If RHS is not the load vector, what is it? If RHS is not the load vector, where is the load vector stored? If RHS is the load vector, why is it changing?

Also while I am here, I don't suppose it might be possible to change the output for the matrix and load vector. Ideally, it would have more significant digits. Also, currently a lot of 0's are printed out though that might just be small numbers that look like zeros because only a few significant digits are printed out.

Typical mass matrix/load vector export:

SELECT
print(matrix)
print(grid)
reorder off
VARIABLES
u
EQUATIONS
div(grad(u))=0
BOUNDARIES
REGION 1
START(0,0) value(u) = 1 LINE TO (1,0) TO (1,1) TO (0,1) value(u) = 0 line TO CLOSE
MONITORS
PLOTS
CONTOUR(u)
END

Alright, have a good day.

Jared

[moderator]

With linear systems, FlexPDE iterates on residuals. The CG solution uses a scaled and preconditioned (actually, pre-and post-), so that minimization in the space of the CG solution is not necessarily minimization in physical space. So after solving the system, the new solution estimate is plugged in to the PDE and a new residual is calculated. Solving against this load produces a new perturbation in the solution estimate. This is repeated until the error tolerance is met or five iterations have been performed (a selectable count, define LINUPDATE). The matrix that is printed is the scaled but not preconditioned version, and it may be rescaled after the first pass. Unfortunately, all these prints were written for debugging purposes, and not for export/import of matrices and solutions.

I have added a facility to control the precision of the matrix output. It will be available in the next maintenance release.

[Jared Barber]

Hello,

Are there any updates on this (i.e. increased precision)?

Also, is there such a facility for this in FlexPDE7? I tried print(matrix), print(grid), reorder off but can't find where FlexPDE7 may have hidden the information.

We are interested in comparing multiple different methods (all employed using FlexPDE) and believe that having a reasonably accurate estimate of the matrix condition number, during all iterations if possible, would help our comparison. Any ideas on how to get that? We could extract it from the matrix printouts but perhaps there is something easier we could try? Perhaps it is already there?

Any info on these things would be helpful. Thanks,

Jared Barber

[moderator]

We have posted developmental version 7.10x3 at http://www.pdesolutions.com/download/xf ... 3win64.exe
In this version, we have done two things:
1) For some reason, the PRINT(MATRIX) selector was never activated in the FlexPDE 7 release version. We have now made it functional. You will get the coupling matrix and RHS every time FlexPDE calls the solver. But reading the output requires some explanation: FlexPDE7 reorders the coupling matrix so that rows are grouped in what we call the I,B,G and K blocks. The I block is the interior nodes, the B block is the (value) boundary nodes, G is the global equations and K is the constraints. The matrix printout is presented as the sixteen matrix blocks II,IB,IG,IK, BI,BB,BG,BK, GI,GB,GG,GK, KI,KB,KG,KK, describing the internal-internal coupling, internal-boundary coupling, etc. There is a reordering map which is not printed, but we can make it available if you need it.
We have added a second command, PRINT(MATRIX15), which will export the same data, but with 15 digits of precision in the values. (This was not the ideal way to do it, but there is no mechanism to specify arbitrary integers in the DEBUG/PRINT handler.)
2) There is a facility which has heretofore been accessible only in our debug versions to compute the eigenvalues of the coupling matrix. We have now made this available in the release version. SELECT DEBUG(MXMODES) to activate this facility. It will print two blocks of data, LambdaR ad LambdaI, the real and imaginary parts of the eigenvalues. It is similar to the DIRECT option, in that it expands the matrix from sparse storage into full storage, so it is a memory hog.

[Jared Barber]

Hello,

Just thought I would note that I am currently running flexpde7 on linux. Hopefully we can get a linux version of 710x3 sometime soon. Thanks,

Jared Barber

[moderator]

We have formally released 7.10 on all platforms. The changes made in 7.10x3 are present in that release.

[Jared Barber]

Oh, ok, thanks. Seems to work. Should be pretty easy to get the condition number using the eigenvalues. Have a good day.

[Jared Barber]

Hello,

I am still having troubles with these facilities. I went ahead and exported the matrix using print(matrix15). I also printed out the eigenvalues using the debug(mxmodes) capability. Finding the eigenvalues did, in fact, take a long time so I read in and looked at the matrix in Matlab. To do some cross-checking I used Matlab to find the condition number of the matrix only to see that the matrix was singular. This contrasts with the estimate we get using the eigenvalues, around 1e5. What I did to construct the matrix was to read in all 16 matrices (II, IB, etc, a lot of which were empty) into sparse 2016x2016 matrices. I then added these matrices together. Looking at their structure, this seemed like the right thing to do.

I have included a zipped file with all that should be needed to produce the matrix and eigenvalue output that I encountered.

I also include below Matlab code that can be used to read in the matrix and eigenvalues. One must change the filename used in "fid = fopen('filename','r');" but everything else should work I think. I tried this out on both linux and windows and got the same results.

Hopefully someone has an idea about this. Thanks and have a good day.

fid = fopen('filename','r');

matrix_names = {'II','IB','IG','IK','BI','BB','BG','BK',...
'GI','GB','GG','GK','KI','KB','KG','KK'};

mystr = fgetl(fid);
while ~strncmpi(' Matrix MATRIX ** length=',mystr,25)
mystr = fgetl(fid);
end
n = sscanf(mystr,' Matrix MATRIX ** length=%d');

for mnc = 1:numel(matrix_names)
curr_matrix_name = matrix_names{mnc};
while ~strncmpi([curr_matrix_name,' matrix'],mystr,9)
mystr = fgetl(fid);
end
mystr = fgetl(fid);
matrix_struc.(curr_matrix_name) = sparse(n,n);
while strncmpi(mystr,'ROW',3)
[row_vals] = sscanf(mystr,'ROW %d OF %d');
i = row_vals(1);
m = row_vals(2); n = row_vals(2);
mystr = fgetl(fid);
while strncmpi(mystr,'(',1)
tmp2 = sscanf(mystr,'(%d) %g ');
inds = tmp2(1:2:end-1);
vals = tmp2(2:2:end);
matrix_struc.(curr_matrix_name)(i,inds) = vals;
mystr = fgetl(fid);
end
end
if ~isfield(matrix_struc,curr_matrix_name)
matrix_struc.(curr_matrix_name) = sparse(m,n);
end

end

BigMatrix = sparse(m,n);
for mnc = 1:numel(matrix_names)
BigMatrix = BigMatrix+matrix_struc.(matrix_names{mnc});
end

% For debugging only
close(figure(1)); figure(1);
for mnc = 1:numel(matrix_names)
subplot(4,4,mnc);
spy(matrix_struc.(matrix_names{mnc}));
end
close(figure(2)); figure(2);
spy(BigMatrix);

%% Now read in the eigenvalues
for ric = 1:2
lc = 0;
while ~contains(mystr,'Lambda')
mystr = fgetl(fid);
end
mystr = fgetl(fid);
k = strfind(mystr,'|');
while ~isempty(k)
mystr = mystr(k+1:end);
tmp = sscanf(mystr,'%g');
nls = numel(tmp);
lambda(lc+[1:nls],ric) = tmp;
lc = lc+nls;
mystr = fgetl(fid);
k = strfind(mystr,'|');
end
end

fclose(fid);