Compute Residuum

Main Forum Old-Forum - do not post here
1

Hi,

After solving my system with the CG method I wanted to compute the relative residuum, but didn’t manage it seems. In my script I do the following

solution = gfu.vec.CreateVector()
solution.data = invpot  * RHS.vec  
gfu.vec.data += solution      

residuum = gfu.vec.CreateVector()
residuum.data = RHS.vec - a.mat*solution

print ("RHS.vec[8:15] = ",RHS.vec[8:15])  
print ("residuum[8:15] = ",residuum[8:15])  
    
norm_residuum = Norm(  residuum)
print ("norm_residuum = ",norm_residuum)
norm_RHS = Norm(RHS.vec)
print ("norm_RHS = ", norm_RHS)        
print ("----------------------- relative residuum = ", norm_residuum/norm_RHS)

And get the following output:

iteration 12 error = 1.9062856897740289e-10
iteration 13 error = 1.529443368104736e-09
iteration 14 error = 4.0000695393358793e-10
iteration 15 error = 4.249320355976443e-11

RHS.vec[8:15] = (6361.74,2.39288e-08)
(13338.1,2.94032e-08)
(13317.2,2.73731e-08)
(6671.39,2.31009e-08)
(-3589.03,-2.1556e-08)
(-15583.3,-3.56244e-08)
(-17718.4,-3.05201e-08)

residuum[8:15] = (1111.11,3.54256e-09)
(2222.22,4.34633e-09)
(2222.22,3.91877e-09)
(1111.11,3.63732e-09)
(-1111.11,-5.11293e-09)
(-2222.22,-5.60155e-09)
(-2222.22,-4.88408e-09)

norm_residuum = 4969.039949981535
norm_RHS = 48245.14687926694
----------------------- relative residuum = 0.1029956435290064

So the CG converged to 1e-12, but my own relative-residual is just indicating 0.1?
Could you please tell me if I made a mistake, and how to avoid that?

One more question: How can I visualize the residuum?

Thanks Jörg

2

Hi Jörg,
If you have dirichlet-dofs you need to project the non-freedofs out, you can do it using the Projector class:

residuum.data = Projector(fes.FreeDofs(), True) * residuum

(it is ok to have the residuum vector on both sides here as the projector works on every entry separately)
To visualize it create a mass matrix and multiply mass times residuum into a GridFunction vector.

Best
christopher

3

Thanks Christopher,

now it works perfectly with this code:
u, v = fesp.TnT()
mass = BilinearForm(fesp)
mass += uvdx
mass.Assemble()

gfu_residuum = GridFunction(fesp)   
gfu_residuum.vec.data = mass.mat * residuum
Draw (gfu_residuum, mesh, "residuum")
4

Sorry, this is the complete code:

solution = gfu.vec.CreateVector()
solution.data = invpot * RHS.vec
gfu.vec.data += solution

residuum = gfu.vec.CreateVector()    
residuum.data = RHS.vec - a.mat*solution
residuum.data = Projector(fesp.FreeDofs(), True) * residuum    

norm_residuum = Norm(  residuum)
print ("norm_residuum = ",norm_residuum)

norm_RHS = Norm(RHS.vec)
print ("norm_RHS = ", norm_RHS)    

print ("----------------------- relative residuum = ", norm_residuum/norm_RHS)

print ("Visualize Residuum")  

#To visualize it create a mass matrix and multiply mass times residuum into a GridFunction vector.
u, v = fesp.TnT()
mass = BilinearForm(fesp)
mass += u*v*dx
mass.Assemble()

gfu_residuum = GridFunction(fesp)   
gfu_residuum.vec.data = mass.mat * residuum
Draw (gfu_residuum, mesh, "residuum")