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