varf/SLEPc matrix definition with arbitrary numbers of equations

General Discussion
1

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’

2

Hi,

my idea would be to use IFMACRO(nGroup,10) ... ENDIFMACRO, and generate the FreeFem code between the IFMACRO-s in some other language (Matlab,Python, etc.) for each number of groups you would like to use. However, I am not an expert, there might be a much better approach.

1 Like
3

This is rather challenging to do in FreeFEM, you can have a look at this paper, Deterministic radiative transfer equation solver on unstructured tetrahedral meshes: Efficient assembly and preconditioning - ScienceDirect, where we solve a problem close to the multi-group diffusion equation. Before turning this into a highly optimized solver, we had a first .edp script that was used to generated a final .edp with any given number of directions (or groups in your case). This is basically what @aszaboa suggested earlier.

1 Like