Hello! I've programmed a structural finite element analysis code in Matlab and have verified it by comparing results with those from the commercial FEA code Abaqus. The Matlab code works fine from a solution accuracy standpoint, but the solve time seems to become unusually large as the size of the system matrix increases, compared to the same analysis in Abaqus.
I've read numerous posts on this forum and a range of other webpages regarding Matlab sparse solvers, but I haven't been able to make any progress on this solution time issue. Maybe I just haven't read the right ones; if so, I apologize in advance for asking a question that has already been answered. Hopefully some people here can help steer me the right way :-) I'll do my best to provide all the relevant information about the Matlab code:
- It uses sparse storage for the system matrix KG and load vector P
- I'm primarily using the \ operator to solve the system: sol = KG \ P
- All analyses are small-displacement and use linear elastic materials. The Abaqus analyses use Abaqus/Standard and are 'linear perturbation' analyses in Abaqus terminology.
- In the small analyses I've run, KG is 7,944 x 7,944. This takes ~4 seconds to solve in Matlab and ~1.2 seconds to solve in Abaqus. Density is 0.0080 via nnz(KG)/prod(size(KG)).
- In the large analyses, KG is 55,479 x 55,479. This takes ~600 seconds (!) to solve in Matlab, but only ~4 seconds to solve in Abaqus. This seems bizarre -- what am I missing? Density is 0.0012.
- For what it's worth, the Abaqus output files state it is using a direct sparse solver (not an iterative method). Here is the key information on this solver as copied from the Abaqus documentation: The direct linear equation solver finds the exact solution to this system of linear equations (up to machine precision). The direct linear equation solver in Abaqus/Standard uses a sparse, direct, Gauss elimination method. The direct sparse solver uses a “multifront” technique that can reduce the computational time to solve the equations dramatically if the equation system has a sparse structure.
- Matlab's iterative methods may offer some potential to speed things up, but from the trials I've done so far, it seems I have to reduce accuracy by an unacceptable amount to get an appreciable reduction in run time.
- I've tried the RCM and AMD reordering functions but have seen no reduction in solution time -- maybe I did this wrong?
abc = symrcm(KG);
testSol.rcm = KG(abc,abc) \ P(abc);
For reference, here are spy plots of the KG matrices for a small model, in the as-built form, and after AMD and RCM. Note that in the as-built cases, there are nonzero terms along the right and bottom 'edges' of KG due to constraint equations. For completeness, the next batch of three spy plots is for a large model (which are fundamentally the same as the small KG matrices). Still scratching my head on what the issue is....
Thanks in advance for any guidance :-)
John
