Building RHS block vector manually

General Discussion
1

Hello all,

I have the following weak formulations:

varf Lhs([dc, dmu], [v, psi]) =

  int2d(Th)((1/chi^2)\*(tau\*eps\*((dx(dmu))\*(dx(v)) + (dy(dmu))\*(dy(v)))) + dc\*v )

  + int2d(Th)( dmu\*psi - (chi^2)\*(3./eps)\*(c\*c)\*dc\*psi - (chi^2)\*eps\*( (dx(dc))\*(dx(psi)) + (dy(dc))\*(dy(psi)) ) );

varf Rhs([dc, dmu], [v, psi]) =

  int2d(Th)((1/chi)\*(tau\*eps\*((dx(mu))\*(dx(v)) + (dy(mu))\*(dy(v))) + (c - c0)\*v) )

  + int2d(Th)( chi\*(mu\*psi - (1/eps)\*(c\*c\*c - c0)\*psi - eps\*((dx(c))\*(dx(psi)) + (dy(c))\*(dy(psi)))) );

And I wish to build the RHS vector manually.

The LHS seems to work, I have built it as follows:

varf vK([u1,u2],[v1,v2]) = int2d(Th) ( dx(u1)*dx(v1) + dy(u1)*dy(v1) );

varf vM([u1,u2],[v1,v2]) = int2d(Th)( u1*v1 );

matrix K = vK(Vh,Vh);

matrix M = vM(Vh,Vh);

varf j(u,v) = int2d(Th)( 3.0*(c*c) * u * v );

matrix J = j(Vh, Vh);

matrix A11 = (1/chi^2)*(tau*eps)*K;

matrix A12 = M;

matrix A21 = M;

matrix A22 = -(chi^2)*J -(chi^2)*K;

matrix A = [

    \[A11, A12\],

    \[A21, A22\]

  \];

And I wish to build the RHS the same way. So far I have:

varf rhsF1([dmu, dc], [v, psi]) = int2d(Th)( (c - c0)*v );

varf rhsF2([dmu, dc], [v, psi]) = int2d(Th)( (tau*eps*( (dx(mu))*(dx(v)) + (dy(mu))*(dy(v)) ) ));

varf rhsG1([dmu, dc], [v, psi]) = int2d(Th)(mu*psi);

varf rhsG2([dmu, dc], [v, psi]) = int2d(Th)(-(1/eps)*(c*c*c - c0)*psi);

varf rhsG3([dmu, dc], [v, psi]) = int2d(Th)( - eps*((dx(c))*(dx(psi)) + (dy(c))*(dy(psi))));

matrix F1 = rhsF1(ME,ME);

matrix F2 = rhsF2(ME,ME);

matrix G1 = rhsG1(ME,ME);

matrix G2 = rhsG2(ME,ME);

matrix G3 = rhsG3(ME,ME);

matrix F = (1/chi)*F1 + (1/chi)*F2;

matrix G = chi*G1 + chi*G2 + chi*G3;

matrix b = [

  [F],

  [G]

];

However, when I try to solve sol = A^-1 * b, I get the following error:

real[int] sol = A^-1 * b; error operator * <18Matrice_Creuse_invIdE>, <14Matrice_CreuseIdE>

Is there an issue with building b like a block matrix? Is there a way to do this with the F and G components as I have them defined?

Thanks in advance.

2

Use a real[int] for the RHS, not a matrix. There are many examples either here or in FreeFEM sources.

3

Is it possible to use it with the manually constructed components like F1,F2, etc.? My goal is to explicitly build the vectors from the variational formulations, without using real[int] b = Rhs(0, ME).

4

Yes, again, that’s how it’s done in many examples.

5

You have to use the syntax real[int] F1 = rhsF1(0, ME) for each of your varfs.
Indeed a(ME,ME) where a is a varf means that you take the bilinear part in a, giving a matrix,
whereas a(0,ME) means that you take the linear part in a, giving a vector.
At the end you build real[int] b=[F,G] and real[int] sol=A^-1*b

1 Like
6

Thank you! That helped a lot, I can now get through the constructions without any errors. However, when I set real[int] sol=A^-1*b; and run the code I get the following error:

ATERROR 1 : VirtualMatrix:: no solver ???

current line = 117

Exec error : MATERROR

– number :1

Exec error : MATERROR

– number :1

err code 8 ,  mpirank

Do I need to manually set the solver? If so, do I set the solver for A,b, or the components F,G ?

I tried set(A, solver = sparsesolver);however, I still receive an error setting up the norm calculations. line 129: real[int] Asol = A*sol;

  current line = 129
Assertion fail : (B->ChecknbColumn(y.N()))
	line :270, in file ./../femlib/RNM.hpp
Assertion fail : (B->ChecknbColumn(y.N()))
	line :270, in file ./../femlib/RNM.hpp
 err code 6 ,  mpirank 0

Thank you for the help, I really appreciate it.

7

You are right, you need to set set(A,solver=sparsesolver);
About your error, try it in two lines:

real[int] Asol(sol.n);
Asol=A*sol;

But maybe it is not that, I can’t say without your full code.

1 Like
8

Unfortunately, that didn’t work. Here is my full code, please take your time, I really do appreciate it.

https://github.com/gart1542/Cahn-Hilliard-FreeFEM-.git

9

This links goes to nothing viewable by me “404 page not found”