# How to impose Robin boundary conditions for a eigenproblem

Hi everyone!
I need to solve the diffusion neutron equaiton in a 2D geometry; In particular i want to solve the eigenproblem. Here below the mathematical description of the problem is shown:

As can be seen from the formulation, the energy is discretized in “G” groups, and the parameter “K_eff” is the eigenvalue.

At this point i would like to impose an “incoming current” = 0 and the icoming current is defined as:

And the vector “J_g” is defined as:

so the final boundary condition should be:

How can i impose this kind of boundary condition for solving the eigenproblem by using the “varf” formulation?

Here is reported the variational form of the problem in my script (NOTICE: here zero FLUX condition is imposed, so a Dirichlet BC, but it is not correct)

``````
func Pk = [P1,P1,P1];
fespace Vh(coremsh,P1);
fespace VhV(coremsh,Pk);
Vh  phi1, phi2, phi3, v1, v2, v3, fisspow, phiTOT;
macro div(D,phi,v) D*(dx(phi)*dx(v) + dy(phi)*dy(v)) // end of macro

varf righthandside([phi1,phi2,phi3],[v1,v2,v3])
= int2d(coremsh)  (CHI[0]*NSF[0]*phi1*v1+CHI[0]*NSF[1]*phi2*v1+CHI[0]*NSF[2]*phi3*v1+CHI[1]*NSF[0]*phi1*v2
+CHI[1]*NSF[1]*phi2*v2+CHI[1]*NSF[2]*phi3*v2+CHI[2]*NSF[0]*phi1*v3+CHI[2]*NSF[1]*phi2*v3+CHI[2]*NSF[2]*phi3*v3);

varf lefthandside([phi1,phi2,phi3],[v1,v2,v3])
= int2d(coremsh)  (div(DFC[0],phi1,v1)+RXS[0]*phi1*v1+div(DFC[1],phi2,v2)+RXS[1]*phi2*v2+div(DFC[2],phi3,v3)+RXS[2]*phi3*v3
- (S1[0]*v1*phi1+S2[0]*v1*phi2+S3[0]*v1*phi3 + S1[1]*v2*phi1+S2[1]*v2*phi2+S3[1]*v1*phi3 + S1[2]*v3*phi1+S2[2]*v3*phi2+S3[2]*v3*phi3))
+ on(0, phi1=0,phi2=0,phi3=0);

real sigma=0;
int nev=1;
real[int] ev(nev); // to store nev eigenvalues
//Vh[int] Vev(nev);
VhV[int] [e1,e2,e3](nev);   // to store nev eigenvectors

matrix A = lefthandside(VhV,VhV, tgv = tgv,solver=GMRES);
matrix B = righthandside(VhV,VhV, tgv = tgv,solver=GMRES);

int kk = EigenValue(A,B, sym=false, sigma=0, value=ev, vector=e1, maxit=0, ncv=0,tol=1e-20);

``````

I think you need to use the `int1d` command as in the thermal conduction example in the documentation.

1 Like

Yes, it is the right command to use for Robin’s BC; It finally turned out that the problem i was facing was about the syntax.

The BC’s expression was something like that :

``````varf problemexpression  = int2d(Th) (...) +int1d(Th,label_i)(1./2.*phi*w);
``````

But it did not work. I finally undestood that, for some reasons, FF seems to reject the usage of “1./2.” directilly inside the BC’s expression, so it is necessary to previously define a real variable, like:

``````real coeff=1./2.;
varf problemexpression  = int2d(Th) (...) +int1d(Th,label_i)(coeff*phi*w);

``````

remark, the eigen value problem must be defined on vector space so your Robin condition must homogeneous.

un exemple:

``````verbosity=1;
mesh Th=square(20,20,[pi*y,pi*x]);
fespace Vh(Th,P2);
Vh u1,u2;

real sigma = 00;  // value of the shift

varf  a(u1,u2)= int2d(Th)(  dx(u1)*dx(u2) + dy(u1)*dy(u2) - sigma* u1*u2 )
+ int1d(Th,1)(u1*u2)
+  on(2,3,4,u1=0) ;  // Boundary condition

varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ; // no  Boundary condition

matrix A= a(Vh,Vh,solver=sparsesolver);
matrix B= b(Vh,Vh,solver=CG,eps=1e-20);

// important remark:
// the boundary condition is make with exact penalisation:
//     we put 1e30=tgv  on the diagonal term of the lock degre of freedom.
//  So take dirichlet boundary condition just on \$a\$ variationnal form
// and not on  \$b\$ variationnanl form.
// because we solve
//  \$\$ w=A^-1*B*v \$\$

int nev=20;  // number of computed eigen valeu close to sigma

real[int] ev(nev); // to store nev eigein value
Vh[int] eV(nev);   // to store nev eigen vector

int k=EigenValue(A,B,sym=true,sigma=sigma,value=ev,vector=eV,tol=1e-10,maxit=0,ncv=0,which="LM"); // which="LM" default
//   tol= the tolerace
//   maxit= the maximal iteration see arpack doc.
//   ncv   see arpack doc.
//  the return value is number of converged eigen value.
``````
1 Like