How to impose Robin boundary conditions for a eigenproblem

General Discussion
1

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:

Screenshot_20230122_203515

Screenshot_20230122_203525
Screenshot_20230122_203525762×141 11.4 KB

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:

image

And the vector “J_g” is defined as:
Screenshot_20230122_205303

so the final boundary condition should be:

Screenshot_20230122_210202

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);
2

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

3

Thanks for your answer!
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);
4

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.