C++

Forum

 
Dirichlet boundary condition in CG

Cplusc (457)
I need to solve a system of linear equation for finite difference method. Actually the problem is to solve the Poisson equation -(d2u/dx2+d2u/dy2)=f(x,y) in which the f=0.25(x^2+y^2) and the domain [-1,1]*[-1,1]. Before applying Dirichlet boundary condition the coefficient matrix is as follows

-1.0 0.5 0.0 0.0 0.5 0.0 0.0 0.0 0.0 
-2.0 0.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.5 -1.0 0.0 0.0 0.50 0.0 0.0 0.0
0.5 0.0 0.0 -2.0 1.0 0.0 0.5 0.0 0.0 
0.0 1.0 0.0 1.0 -4.0 1.0 0.0 1.0 0.0
0.0 0.0 0.5 0.0 1.0 -2.0 0.0 0.0 0.5
0.0 0.0 0.0 0.5 0.0 0.0 -1.0 0.5 0.0 
0.0 0.0 0.0 0.0 1.0 0.0 0.5 -2.0 0.5 
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 -1.0


and the right hand side is 
{0.33,0.5,0.16,0.5,1,0.5,0.16,0.5,0.33}


nodes {1,2,3,4,6,7,8,9} are the boundary points and corresponding boundary value for them are {0.5,0.25,0.5,0.25,0.25,0.5,0.25,0.5}.
I applied the boundary condition to both coefficient matrix and right hand side as follows

1 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 
0 1 0 1 -4 1 0 1 0
0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 1


rhs={0.5,0.25,0.5,0.25,1,0.25,0.5,0.25,0.5}


the I used the conjugate gradient solver for that.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
 
void Blas::CG(VecVecDbl_t& A, VecDbl_t& x, VecDbl_t& b, size_t iter_max) {

	size_t N = A.size();
	VecDbl_t r(N, 0.0);
	VecDbl_t rn(N, 0.0);
	VecDbl_t p(N, 0.0);
	VecDbl_t pn(N, 0.0);
	VecDbl_t Ap(N);
	x.resize(N);
	VecDbl_t xn(N);
	double tol = 1e-6;

	for (int i = 0; i < b.size(); i++)
	{
		 r[i] = p[i] = b[i];
				 
	}
	size_t iter = 0;
	while (iter < 1)
	{
		double rtr = 0.0;
		double ptAp = 0.0;
		rtr = dotproduct(r, r);
		MPI_MxV(A, p, Ap);
		ptAp = dotproduct(p, Ap);

		double alpha = rtr / (ptAp);
		
		for (int i = 0; i < N; i++)
		{
			x[i] += alpha * p[i];
			rn[i] = r[i]-alpha*Ap[i];

		}
		double rntrn = 0.0;
		rntrn = dotproduct(rn, rn);

		if (vectorNorm(r) < tol) break;
		double beta = rntrn / rtr;

		for (int i = 0; i < N; i++)
		{
			pn[i] = rn[i] + beta * p[i];
		}
		iter++;

		p = pn;
		r = rn;
	}

}

Then I don't understand why the boundary condition are not set at the first iteration. can anyone tell me where am I wrong?
Last edited on
TheIdeasMan (6817)
Please don't start new topics on the same subject, asking the exact same question over again. What was wrong with the answer in your previous topic?
Cplusc (457)
Without knowing what physical system you are trying to solve for it's impossible to comment further.

I tried to explain more about the physical system. and also I initialized the boundary condition before iteration loop, but I have the same problem, the dirichlet boundary value are not set in the solution vector in the first iteration.
TheIdeasMan (6817)
Yes, but give people time to answer, don't start a new topic. Otherwise it will like all your other posts, we go round & round in circles.
Cplusc (457)
ok,sure
TheIdeasMan (6817)
I don't know anything about the math for this problem, just giving some advice for you to help yourself.

resabzr wrote:
.... the dirichlet boundary value are not set in the solution vector in the first iteration.


Instead of just saying values are not set, provide some evidence of it. Get your program to print the values on the first iteration. Post them here. Saying the values are not set doesn't seem to make sense: lastchance was saying that the Dirichlet boundary condition values are set at the beginning and remain constant, how and why is it different for your case?

May be it's only a perception that the values are wrong?

Have you used a debugger? It's an ideal tool for when values are wrong, or there is a logic problem.

In general, about posting:

Imagine every time you post, it's like asking your boss for help. Is your boss going to be happy that you have asked 160 times for help? What can you do to help yourself? Here are some tips:

If using an external library, read the documentation for it;
If the program doesn't compile, look at cppreference.com for help;
If the program does compile, and there is a value or logic problem, use a debugger;
Look at wiki pages about the subject;
Search the web in general, someone is bound to have the same problem, look at their code;
If you do post on a forum, make it as easy possible for those who are going to help you, provide code that can be compiled, provide input and expected output;
Show that you have implemented or at least tried suggestions from the forum.

Are you aware that the while loop in the code above only runs once?

Does vectorNorm(r) return a double? Maybe it is poorly named, returns the magnitude?
lastchance (6980)
@resabzr,
You appear to be trying to set fixed boundary values (which is all that is meant by Dirichlet boundary conditions) by using your matrix equation. i.e. you have lines in your matrix with a single 1 on the diagonal so that you set that value to the equivalent value on the RHS:
( . . . . . )(.)   (.)
( 0 1 0 0 0 )(x)   (c)
( . . . . . )(.)   (.)
( . . . . . )(.) = (.)
( . . . . . )(.)   (.)
( . . . . . )(.)   (.)

will set x to c when you have finished inverting the matrix.

The fact that you can write the problem down mathematically this way ... does not mean you should code it like that. Inverting matrices is a very expensive operation - you should not increase the size of your matrix by including trivial operations (like setting fixed values), because it will make the matrix inversion unnecessarily more costly.

You claim to be using an iterative algorithm. That means you require a number of iterations to invert your matrix system. In a really badly-preconditioned system it could take N iterations for an NxN matrix (although that is unlikely for methods like conjugate gradients or GMRES). There is absolutely no reason that any particular element in your solution vector will reach its final value after only one iteration. The only solution procedures that would do that here are methods like Gauss-Seidel, and they won't match conjugate gradients for large matrices.

The correct place for fixed boundary conditions in a sparse matrix system is on the RHS ... because those values are already known. On the left, in the solution vector, you should only have points whose values are yet to be found. That way, your matrix is as small as possible and the matrix inversion cost is minimal.
Last edited on
Cplusc (457)
Thank you all for answers.
Last edited on
oggin (51)

Row echelon will do this with full pivoting:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
 
 


#include<iostream>  
#include<vector> 
#include<cmath> 
#include <iomanip>

using namespace std;  

inline void pivot(vector<vector<double>> &b,vector<double>&r,int num,int n)
{
	int p1,p2,g;
	for (p1 = num;p1<=n-1;p1++) 
		{
				for (p2 = p1+1;p2<=n;p2++)
			{
				if (  abs(b[p1][num]) < abs(b[p2][num])  )
			    {
			swap(r[p1],r[p2]);
			for(g=0;g<=n;g++)
			{
				swap (b[p1][g],b[p2][g]);
			}
	    	}	
	    }
	}	
}

void gauss(vector<vector<double>> b,vector<double> r, vector<double> &ans)
{
	int n=b.size(),p1,p2,num,k,row,g,z,j;
	double f;
	ans.resize(n);
	n=n-1;
for(k=0;k<=n-1;k++)
{
	pivot(b,r,k,n);
	for(row=k;row<=n-1;row++)
	{
		 if (b[row+1][k]==0) {break;};
		  f=b[k][k]/b[row+1][k];
		  r[row+1]=r[row+1]*f-r[k];
		 for(g=0;g<=n;g++)
		 {
		 	 b[row+1][g]=b[row+1][g]*f-b[k][g];
	     }
	}	
}

//back substitute
for(z=n;z>=0;z--)
{
	 ans[z]=r[z]/b[z][z];
	 for(j=n;j>=z+1;j--)
	 {
	 	 ans[z]=ans[z]-(b[z][j]*ans[j]/b[z][z]);
	 }
}
}//gauss


void matmultvect(vector<vector<double>> m1,vector<double> m2, vector<double> &ans)
{
	int rows=m1.size(),r,k;
	double rxc;
	ans.resize(rows);
	rows=rows-1;
	for(r=0;r<=rows;r++)
	{
		rxc=0;
		for(k=0;k<=rows;k++)
		{
			rxc=rxc+m1[r][k]*m2[k];
		}
	ans[r]=rxc;	
	}	
}


void print(vector<double> v)
{
	int n=v.size(),i;
	string sp;
	for(i=0;i<n;i++)
	{
		if(v[i]>=0) sp="+"; else sp="";
		cout<<setprecision(15)<<sp<<v[i]<<endl;
	}	
}


int main()
{
	
	int n;
vector<vector<double>> d 
{
{-1.0,0.5,0.0,0.0,0.5,0.0,0.0,0.0,0.0},
{-2.0,0.5,1.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.5,-1.0,0.0,0.0,0.50,0.0,0.0,0.0},
{0.5,0.0,0.0,-2.0,1.0,0.0,0.5,0.0,0.0},
{0.0,1.0,0.0,1.0,-4.0,1.0,0.0,1.0,0.0},
{0.0,0.0,0.5,0.0,1.0,-2.0,0.0,0.0,0.5},
{0.0,0.0,0.0,0.5,0.0,0.0,-1.0,0.5,0.0},
{0.0,0.0,0.0,0.0,1.0,0.0,0.5,-2.0,0.5},
{0.0,0.0,0.0,0.0,0.0,0.5,0.0,0.5,-1.0},
};


    
 vector<double> r={0.5,0.25,0.5,0.25,1,0.25,0.5,0.25,0.5};


vector<double> ret;

gauss(d,r,ret);

cout<<"solution:"<<endl;
print(ret);

cout <<"  "<<endl;
// multiply back to check
vector<double>ret2;
matmultvect(d,ret,ret2);
cout<<"check, should be {0.5,0.25,0.5,0.25,1,0.25,0.5,0.25,0.5:"<<endl;
print(ret2);

cout <<"  "<<endl;
cout <<"Press return to end . . ."<<endl; 
cin.get();	  
	
}
Cplusc (457)
Thanks
Topic archived. No new replies allowed.