Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.
Equation in matrix form for use for MATLAB Livelink
Posted 29 apr 2015, 07:07 GMT-4 Interfacing, Parameters, Variables, & Functions, Studies & Solvers Version 5.2a 9 Replies
Please login with a confirmed email address before reporting spam
This seems to be a pretty simple question but I am struggling to find any help on the topic!
I have recently started using COMSOL with the livelink to write a Model Order Reduction in MATLAB. The goal is to use COMSOL to generate the system matrices and eventually solve the problem in MATLAB.
To this end, the command mphmatrix should be of great help. The help page for this function tells us that we can retrieve following system matrices:
'K' Stiffness matrix
'L' Load vector
'M' Constraint vector
'N' Constraint Jacobian
'D' Damping matrix
'E' Mass matrix
'NF' Constraint force Jacobian
'NP' Optimization constraint
Jacobian(*)
'MP' Optimization constraint vector(*)
'MLB' Lower bound constraint vector(*)
'MUB' Upper bound constraint vector(*)
'Kc' Eliminated stiffness matrix
'Lc' Eliminated load vector
'Dc' Eliminated damping matrix
'Ec' Eliminated mass matrix
'Null' Constraint null-space basis
'Nullf' Constraint force null-space
matrix
'ud' Particular solution ud
'uscale' Scale vector
Moreover with the command mphgetu, one can retrieve solution vector U and its first derivative Udot at any given time.
The problem is nowhere can I find what linear system is COMSOL solving. I can only fathom it is:
E*Udotdot + D*Udot + K*U == L
What is even more surprising is that the model has 2315 degrees of freedom (dofs) plus 6 internal dofs. Matrices E, D and K are 2315x2315, vector L is 2315x1. Very well. However U as retrieved with command mphgetu is 2321x1 ! U seems to include the 6 internal dofs.
My question is, where can I find help/documentation on what linear system COMSOL solves, and where can one find what is the correlation between U and dofs. mphxmeshinfo seems promising but the documentation for that function is vague.
Many many thanks for all help.
Best regards,
Guillaume
Hello Guillaume Tyson
Your Discussion has gone 30 days without a reply. If you still need help with COMSOL and have an on-subscription license, please visit our Support Center for help.
If you do not hold an on-subscription license, you may find an answer in another Discussion or in the Knowledge Base.
Please login with a confirmed email address before reporting spam
Could you express for me in details!
Thanks a lot!
Please login with a confirmed email address before reporting spam
So if you check in the COMSOL documentation you'll find in "About the Time-Dependent Solver" exactly what equation COMSOL solves.
I'm not quite sure I understand what the "eliminated" matrices are but to answer your question, what I am sure of is that with the matrices from above, the equation is written:
E*Udotdot + D*Udot + K*U == F
Where F is a new vector containing the flux terms for the different dofs. There are two ways to calculate F:
-Either you use the complete system residuals at each time step that are actually in matrix L. This means that at any time step you would calculate F by:
F = E*Udotdot + D*Udot + K*U + L
This is very slow.
-Or you assemble F yourself from the boundary conditions of your problem. This can be tricky and I had to ask support for help but it works very well and is very fast.
Also maybe your problem comes from the presence of additional dofs. You can turn them off by clicking the "Show" button in the Model Builder then ticking "Discretization". You'll see that in your model interface, under the "Equation" drop-down menu a "Discretization" drop-down menu has appeared. "Compute boundary fluxes" should be checked, turn this option off and then the vector U in Matlab should have as many lines as the system matrices.
Hope this helps,
Regards,
Guillaume
Please login with a confirmed email address before reporting spam
the problem i am concerned is the mass or damping matrix.
for example, briefly, how to extract the mass matrix of poisson equation with certain boundary condition.
Could you give me future help ?
Best regards!
- Zhu Shuai
Please login with a confirmed email address before reporting spam
The problem is that you can't really extract the boundary conditions. After the call MAT = mphmatrix(model, 'sol1', 'out', {'L'}), MAT.L will be the load vector. As soon as you specify a solution number ( ie mphmatrix(model, 'sol1', 'out', {'L'}, i) for the i-th solution time) MAT.L will become your equation residuals.
You should probably aim to rebuild the the load vector manually in MATLAB.
Regards,
Guillaume
Please login with a confirmed email address before reporting spam
So, we can also operate the boundary condition for the mass matrix in the same way. Is there any easy way to implement it via Comsol ?
Best Regards!
---ZhuShuai
Please login with a confirmed email address before reporting spam
Have you found the reason of this added 'internal DOFs'? Or did you get the correlation between the solution vector U (via getU) and other vectors? I need to extract the U, modify it and then assign it back to variables, thus I have to know how the U is formatted.
Regards
Hao
Please login with a confirmed email address before reporting spam
You can use the function mphxmeshinfo to get information the matrices returned by mphmatrix and the internal degrees of freedom. The data returned by mphxmeshinfo is documented in the Programming Reference Manual under XMeshInfo in the chapter on solvers.
Rregards
Lars Gregersen
Comsol Denmark
Please login with a confirmed email address before reporting spam
Please login with a confirmed email address before reporting spam
Hi,
The problem is that you can't really extract the boundary conditions. After the call MAT = mphmatrix(model, 'sol1', 'out', {'L'}), MAT.L will be the load vector. As soon as you specify a solution number ( ie mphmatrix(model, 'sol1', 'out', {'L'}, i) for the i-th solution time) MAT.L will become your equation residuals.
You should probably aim to rebuild the the load vector manually in MATLAB.
Regards,
Guillaume