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:

Screenshot_20230122_203515

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

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

1 Like

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