Dear FF++ developers,

I am wondering how I can define a macro or a function that takes in input the number of equations and automatically defines the variational forms and the matrices for the system solution.

More specifically, my application involves the solution of the multi-group diffusion equation for neutrons, which is basically a system of diffusion equations per each energy group used to discretise the energy axis.

In the following, you can find my definition for the two-group case. First, I define the variational forms and the associated matrices for each term of the neutron diffusion equation, then I pass the final matrices to EPSSolve for the solution with SLEPc.

My goal would be to select an arbitrary number of groups (e.g., 42) without having to define my equations by hand, but I cannot see a way to do that (e.g., *varf scattering([phi1,phi2],[v1,v2])* would require *varf scattering([phi1,phi2,…phi42],[v1,v2,…,v42])*, which is absolutely not flexible).

```
// fespace definition
func Pk = [P1, P1];
fespace Vh(coremsh,P1);
fespace VhV(coremsh,Pk);
Vh phi1, phi2, v1, v2, fisspow, phiTOT;
macro div(D,phi,v) D*(dx(phi)*dx(v) + dy(phi)*dy(v)) // end of macro
// define variational forms
// the capitalised variables are the nuclear properties per each energy group, they are simple P0 fields
varf leakage([phi1,phi2],[v1,v2]) = int2d(coremsh) (div(DFC[0], phi1, v1)+div(DFC[1], phi2, v2))
+on(1,2,3,4, phi1=0, phi2=0);
varf removal([phi1,phi2],[v1,v2]) = int2d(coremsh) (RXS[0]*phi1*v1 + RXS[1]*phi2*v2);
varf scattering([phi1,phi2],[v1,v2]) = int2d(coremsh) (S12[0]*phi1*v2);
varf fission([phi1,phi2],[v1,v2]) = int2d(coremsh) (NSF[0]*phi1*v1+NSF[1]*phi2*v1);
// build operators
matrix<real> L = leakage(VhV, VhV, tgv = tgv);
matrix<real> S = scattering(VhV, VhV, tgv = tgv);
matrix<real> R = removal(VhV, VhV, tgv = tgv);
matrix<real> F = fission(VhV, VhV, tgv = tgv);
matrix<real> LHS = L+R-S;
matrix<real> RHS = F;
Mat A(LHS, clean = true);
Mat B(A, RHS, clean = true);
/* ----------------------- SOLVE EIGENVALUE PROBLEM ------------------------------------------------------- */
real sigma=0;
int nev=1;
real[int] ev(nev); // to store nev eigenvalues
VhV<real>[int] [e1,e2](nev); // to store nev eigenvectors
string ssparams = // Parameters for coremshe distributed EigenValue solver
" -eps_nev " + nev + // Number of eigenvalues
" -eps_type krylovschur" +
" -eps_target "+ sigma + // Shift value
" -st_type sinvert " +
" -st_pc_type lu " +
" -st_pc_factor_mat_solver_type mumps" +
" -eps_view"
;
int k = EPSSolve(A, // matrix OP = A − sigma*B
B, //
vectors = e1, // Array to store coremshe FEM-EigenFunctions
values = ev, // Array to store coremshe EigenValues
sparams = ssparams // Parameters for coremshe distributed EigenValue solver
);
```

Can I use a for loop to define an array of varf or something similar? How would you proceed? It is very difficult for me to imagine a way to do that automatically with the current syntax…

Thank you so much for your support and help!

Nicolo’