esys-escript

residual unbalanced force and projector

Asked by ceguo

Hi,

I have two questions wishing anyone could help.

One is how to calculate the residual unbalanced force which is used to judge convergence when nonlinear problems are solved. Usually it's calculated from the nodal force in weak form. Then how to calculate it when strong form is used as in Escript?

The second question is about the use of Projector. Here I use it to project the deviatoric stress from the gaussian points to the nodes for visualization. I hope Projector does serve this purpose. Then how are the quantities projected from integration points to nodal points?

Thanks in advance.

Ning

Question information

Language:
English Edit question
Status:
Solved
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Solved by:
Lutz Gross
Solved:
Last query:
Last reply:
Revision history for this message
Lutz Gross (l-gross) said : 		#1

ad Q1: Escript does not use the strong form. You can use the getResidual method to get the residual forces at each node.

ad Q2: The projector solves the PDE "D*u=Y" (in escript notation) where D=1, Y=quantity on integration points and - as return - u is the quantity on points. By default lumping is used on the stiffness matrix. This means for
point i one gets

u_i = sum_e ( sum_j Y[j,e]*w[j,e]*phi[i,j,e]) / sum_e (sum_j w[j,e]*phi[i,j,e])

where Y[j,e] is the value of Y at the j-th integration point in element e, w[j,e] is the integration weight,
 at the j-th integration point in element e and phi[i,j,e] is the shape function for point i at at the j-th integration point in element e.

Revision history for this message
ceguo (hhh-guo) said : 		#2

Thanks for the answer.

But I still have problem on calculating the residual at each point. When I call mypde.getResidual() I encountered a RuntimeError: FinleyAdapterException: Paso_SystemMatrix_MatrixVector: balanced matrix is not supported. (mypde is constructed from mypde=LinearPDE(mydomain))

Also I don't think this will return what I want. My problem is a nonlinear problem which I wrote a Newton-Ralphson scheme to solve. After solving

-(A_{ijkl} u_{k,l})_{,j} = -X_{ij,j} + Y_i, where A_{ijkl} is the tangent operator, -X_{ij} is the current stress and Y_i is the body force.

I applied grad(u) to the constitutive model and got true stress. And now I want to get the residual to see if the problem is converged and if any further iterations are needed. It is something like R_i=-X'_{ij,j}+Y_i, where -X'_{ij} is updated stress.

So my question is how to calculate R_i and extrapolate it to the residual nodal force. Thanks.

Ning

Revision history for this message
Best Lutz Gross (l-gross) said : 		#3

Ops! Once you have solved the PDE the stiffness matrix is modified which stops you to calculate the residual. I made a note to fix this.

To calculate some sort of residual of the current approximation you can use the LinearPDE.getRightHandSide():
If you run a Newton-Ralphson scheme by solving

-(A_{ijkl} u_{k,l})_{,j} = -X_{ij,j} + Y_i,

to get the correction u then you can use

    mypde.setValue(A=..., X=..., Y=...)
    r=mypde.getRightHandSide()
    norm_r=Lsup(r)

r is in fact given on the nodes of the FEM mesh.

If phi_k is the shape function for node k then r_ik = int_{domain} ( phi_k,j X_{ij) + phi_k Y_i) dx.
As phi_k has compact support (and ignoring surface forces) it is in fact
r_ik = int_{domain} phi_k (Y_i- X_{ij,i} ) dx which is your residual force (you may want to rescale this by (int_{domain} phi_k)^{-1} )

Revision history for this message
ceguo (hhh-guo) said : 		#4

Thanks.
The explanation is clear. I just want to make sure that in my problem, the convergence criteria for the Newton-Raphson scheme can be set as:

norm_r=Lsup(r) < tolerance x norm_r_ref,

where norm_r_ref can be the initial residual norm.

This is reasonable, right?

Ning

Revision history for this message
Lutz Gross (l-gross) said : 		#5

If the elements in the mesh have roughly the same size, this is reasonable as int_{domain} phi_k dx ~ const.
If there is extreme local refinement, it would be better to use Lsup(r_k/int_{domain} phi_k dx).
You can get r_k/int_{domain} phi_k dx if you need to (e.g. by using a second instance of the PDE class with Lumping). It is on the to do list that the LinearPDE class can calculate this.

Revision history for this message
ceguo (hhh-guo) said : 		#6

Thanks Lutz Gross, that solved my question.