H1 error for vector function problem

[Hemanta822]

I was trying to compute the H1 error for a vector problem in 2D. As my exact solution, I have a 2D vector u_ex = [u1, u2]. I created a 2*2 gradient matrix and used MatrixFunctionCoefficient to get the coefficients. But when I try to compute the H1 error using u.ComputeH1Error(exact_u_coeff, exact_grad_coeff), I cannot do so because ComputeH1Error() takes gradient from VectorFunctionCoefficient. I would appreciate it if someone could help me or give me some direction to deal with this.

[vladotomov]

Hi @Hemanta822,

I think ComputeH1Error(Coefficient *, VectorCoefficient *, ...) assumes scalar fields. The vector version is not implemented. If you want to implement and add ComputeH1Error(VectorCoefficient *exsol, MatrixCoefficient *exgrad, ...) to MFEM it would be great.

With the existing methods you can split your vector function into scalar GridFunctions and compute H1Error for each.

[Hemanta822]

Hi @vladotomov, I am working on example 5 (with different function space). And I am simplifying my problem by setting both components of displacement (u) same (and alpha = 0, so that I do not see any impact of pressure). I expected to see my x = [x0, x1] such that x0 = x1. However, my x0 and x1 have many same numbers but in the different positions (for example in the following mentioned solution, 0.707232 and 0.70718 have swapped their positions in x0 and x1). I also tried to play with the ordering in the function space as:

//FiniteElementSpace *R_space = new FiniteElementSpace(mesh, fec, dim, Ordering::byVDIM);
   FiniteElementSpace *R_space = new FiniteElementSpace(mesh, fec, dim, Ordering::byNODES); 

I got my displacement x as follows:

Solution (x):
0 0 0 0 0 0 0.500526 0.707232
0.499646 8.65869e-17 0 0.70718 1.00017 0.70718 1.22452e-16 0
0.499646 0.707232 0.500526 8.65869e-17 0 8.65869e-17 1.22452e-16 8.65869e-17
1.49961e-32 0 0.270211 0.270813 0 0 0.382735 0.653444
0.923834 0.923895 0.653088 0.38267 0 0.653361 0.270385 0
0.271093 4.68605e-17 0 0.653218 0.382735 1.13131e-16 0.923834 0.923895
0.653088 0.653444 0.270211 1.13131e-16 0.270813 4.68605e-17 4.68605e-17 0.38267
1.13131e-16 0.270385 1.13131e-16 0.653361 0.653218 0 0.271093 0
4.68605e-17 0.145539 0.353019 0.853213 0.35307 0.353767 0.147396 0.35394
0.853516 0.853213 0.353019 0.145539 0.35307 0.353767 0.853516 0.35394
0.147396 0 0 0 0 0 0 0.500526
0.70718 0.499646 8.65869e-17 0 0.707232 1.00017 0.707232 1.22452e-16
0 0.499646 0.70718 0.500526 8.65869e-17 0 8.65869e-17 1.22452e-16
8.65869e-17 1.49961e-32 0 0.270813 0.270211 0 0 0.38267
0.653088 0.923895 0.923834 0.653444 0.382735 0 0.653218 0.271093
0 0.270385 4.68605e-17 0 0.653361 0.38267 1.13131e-16 0.923895
0.923834 0.653444 0.653088 0.270813 1.13131e-16 0.270211 4.68605e-17 4.68605e-17
0.382735 1.13131e-16 0.271093 1.13131e-16 0.653218 0.653361 0 0.270385
0 4.68605e-17 0.145539 0.35307 0.853213 0.353019 0.35394 0.147396
0.353767 0.853516 0.853213 0.35307 0.145539 0.353019 0.35394 0.853516
0.353767 0.147396

Solution0 (x0):
0 0 0 0 0 0 0.500526 0.707232
0.499646 8.65869e-17 0 0.70718 1.00017 0.70718 1.22452e-16 0
0.499646 0.707232 0.500526 8.65869e-17 0 8.65869e-17 1.22452e-16 8.65869e-17
1.49961e-32 0 0.270211 0.270813 0 0 0.382735 0.653444
0.923834 0.923895 0.653088 0.38267 0 0.653361 0.270385 0
0.271093 4.68605e-17 0 0.653218 0.382735 1.13131e-16 0.923834 0.923895
0.653088 0.653444 0.270211 1.13131e-16 0.270813 4.68605e-17 4.68605e-17 0.38267
1.13131e-16 0.270385 1.13131e-16 0.653361 0.653218 0 0.271093 0
4.68605e-17 0.145539 0.353019 0.853213 0.35307 0.353767 0.147396 0.35394
0.853516 0.853213 0.353019 0.145539 0.35307 0.353767 0.853516 0.35394
0.147396

Solution1 (x1):
0 0 0 0 0 0 0.500526 0.70718
0.499646 8.65869e-17 0 0.707232 1.00017 0.707232 1.22452e-16 0
0.499646 0.70718 0.500526 8.65869e-17 0 8.65869e-17 1.22452e-16 8.65869e-17
1.49961e-32 0 0.270813 0.270211 0 0 0.38267 0.653088
0.923895 0.923834 0.653444 0.382735 0 0.653218 0.271093 0
0.270385 4.68605e-17 0 0.653361 0.38267 1.13131e-16 0.923895 0.923834
0.653444 0.653088 0.270813 1.13131e-16 0.270211 4.68605e-17 4.68605e-17 0.382735
1.13131e-16 0.271093 1.13131e-16 0.653218 0.653361 0 0.270385 0
4.68605e-17 0.145539 0.35307 0.853213 0.353019 0.35394 0.147396 0.353767
0.853516 0.853213 0.35307 0.145539 0.353019 0.35394 0.853516 0.353767
0.147396

Please let me know how I can fix it and if there is any easier way to do it. Thank you so much!

[vladotomov]

I can't say from the given info. What is alpha, I don't see such variable in ex5.cpp?
If only few of the numbers are different, my first guess would be boundary conditions.

[Hemanta822]

Thank you for your response @vladotomov. In example 5, if I want to compute H1 error for u, how would I do it? After we solved using MINRES, we used our solution x to get grid function u as:

u.MakeRef(R_space, x.GetBlock(0), 0);

This includes both the components of u. From your first comment, I thought I could just split x into x0 and x1 and get gridfunctions u1 and u2. But, I think I am not getting correct x0 and x1. The boundary conditions are correctly imposed. My all other errors where I didn't split x and used gridfunction u, I am getting correct optimal rates. So, there must be something wrong with the way I split x.

I did the following:

int x_length = x.GetBlock(0).Size();
int x_half = x_length/2;
 
for (int i = 0; i < x_half ; i++)
    {
         x0[i]=x[i];
    }
 
   for (int i = 0; i < x_half ; i++)
    {
         x1[i]=x[i+x_half];
    }

Please let me know your thought. Thank you!

[vladotomov]

Are you solving ex5 with the RT space? The above won't extract the components of an RT function correctly, that's not how the vector unknowns are arranged in the solution vectors for RT. See this comment for a way to project the components to an L2 vector space. For L2, the above extraction should work.

[Hemanta822]

Hi @vladotomov, I am using H1_FECollection,

FiniteElementCollection *fec = new H1_FECollection(order + 1, dim);
FiniteElementSpace *R_space = new FiniteElementSpace(mesh, fec, dim, Ordering::byVDIM); 

Since I am having dropped H1 rates while using these splitted x0 and x1. I checked with my L2 rates.

I am getting optimal rate 3.00, when I compute the L2 norm using the existing functionality in MFEM WITHOUT SPLITTING x as follows:

u.MakeRef(R_space, x.GetBlock(0), 0);
real_t l2_err_u  = u.ComputeL2Error(ucoeff, irs);

However I see high drop in my L2 rate (similar to what I observe in H1 rate), if I split my x as mentioned above and compute L2 error using components as follows:

double l2_err_0 = u0.ComputeL2Error(ucoeff0, irs);
double l2_err_1 = u1.ComputeL2Error(ucoeff1, irs);
   
double l2_err_u = sqrt(l2_err_0*l2_err_0 + l2_err_1*l2_err_1);

where u1 and u2 are gridfunctions created using x0 and x1 (x = [x0, x1]).

This is why I am guessing if the way I splitted the solution doesn't work. Please advise. Thank you!

[vladotomov]

I see, interesting, can you send me your version of ex5.cpp and I'll take a look.

[vladotomov]

Hi @Hemanta822,

Make sure you call SetTime for ucoeff0 and ucoeff1 as the exact solutions are time-dependent.

[Hemanta822]

Hi @vladotomov, Thank you so much for helping with this. It gives me the expected rates now. Thank you again!