diff --git a/src/pdes/neutron_sn/bc.c b/src/pdes/neutron_sn/bc.c new file mode 100644 index 00000000..184b10a5 --- /dev/null +++ b/src/pdes/neutron_sn/bc.c @@ -0,0 +1,133 @@ +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: boundary conditions + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#include "feenox.h" +#include "neutron_sn.h" + +int feenox_problem_bc_parse_neutron_sn(bc_data_t *bc_data, const char *lhs, char *rhs) { + + // TODO: should this be the default BC? + if (strcmp(lhs, "vacuum") == 0 || strcmp(lhs, "null") == 0) { + // "null" is supported for compatibility with diffusion + bc_data->type_math = bc_type_math_dirichlet; + bc_data->set_essential = feenox_problem_bc_set_neutron_sn_vacuum; + bc_data->dof = -1; + + } else if (strcmp(lhs, "mirror") == 0 || strcmp(lhs, "symmetry") == 0) { + bc_data->type_math = bc_type_math_dirichlet; + bc_data->dof = -1; + bc_data->set_essential = feenox_problem_bc_set_neutron_sn_mirror; + + // TODO: albedo + // TODO: white current + // TODO: individual scalar fluxes as a function of theta & phi + } else { + feenox_push_error_message("unknown neutron_sn boundary condition '%s'", lhs); + return FEENOX_ERROR; + } + + bc_data->space_dependent = feenox_expression_depends_on_space(bc_data->expr.variables); + bc_data->nonlinear = feenox_expression_depends_on_function(bc_data->expr.functions, feenox.pde.solution[0]); + + if (bc_data->nonlinear && bc_data->type_math == bc_type_math_dirichlet) { + feenox_push_error_message("essential boundary condition '%s' cannot depend on phi", rhs); + return FEENOX_ERROR; + } + + feenox_call(feenox_expression_parse(&bc_data->expr, rhs)); + + // TODO: check consistency, non-linearities, etc. + + return FEENOX_OK; +} + + +// these virtual method fill in the dirichlet indexes and values with bc_data +int feenox_problem_bc_set_neutron_sn_vacuum(bc_data_t *this, element_t *e, size_t j_global) { + +#ifdef HAVE_PETSC + double outward_normal[3]; + feenox_call(feenox_mesh_compute_outward_normal(e, outward_normal)); + for (unsigned n = 0; n < neutron_sn.directions; n++) { + if (feenox_mesh_dot(neutron_sn.Omega[n], outward_normal) < 0) { + // if the direction is inward set it to zero + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + feenox_call(feenox_problem_dirichlet_add(feenox.pde.mesh->node[j_global].index_dof[dof_index(n,g)], 0)); + } + } + } +#endif + + return FEENOX_OK; +} + + +int feenox_problem_bc_set_neutron_sn_mirror(bc_data_t *this, element_t *e, size_t j_global) { + +#ifdef HAVE_PETSC + double outward_normal[3]; + double reflected[3]; + double Omega_dot_outward = 0; + double eps = feenox_var_value(feenox.mesh.vars.eps); + + // TODO: mark the BC as dependent on the normal and compute it in the caller + feenox_call(feenox_mesh_compute_outward_normal(e, outward_normal)); +// printf("element %lu node %lu outward normal %.3f %.3f %.3f\n", e->tag, j_global, outward_normal[0], outward_normal[1], outward_normal[2]); + for (unsigned n = 0; n < neutron_sn.directions; n++) { + if ((Omega_dot_outward = feenox_mesh_dot(neutron_sn.Omega[n], outward_normal)) < 0) { + // if the direction is inward then we have to reflect it + // if Omega is the incident direction with respect to the outward normal then + // reflected = Omega - 2*(Omega dot outward_normal) * outward_normal + + for (int m = 0; m < 3; m++) { + reflected[m] = neutron_sn.Omega[n][m] - 2*Omega_dot_outward * outward_normal[m]; + } + + unsigned int n_prime = 0; + for (n_prime = 0; n_prime < neutron_sn.directions; n_prime++) { + if (fabs(reflected[0] - neutron_sn.Omega[n_prime][0]) < eps && + fabs(reflected[1] - neutron_sn.Omega[n_prime][1]) < eps && + fabs(reflected[2] - neutron_sn.Omega[n_prime][2]) < eps) { + break; + } + } + + if (n_prime == neutron_sn.directions) { + feenox_push_error_message("cannot find a reflected direction for n=%d (out of %d in S%d) for node %d", n, neutron_sn.directions, neutron_sn.N, feenox.pde.mesh->node[j_global].tag); + return FEENOX_ERROR; + } + + double *coefficients = NULL; + feenox_check_alloc(coefficients = calloc(feenox.pde.dofs, sizeof(double))); + + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + coefficients[dof_index(n,g)] = -1; + coefficients[dof_index(n_prime,g)] = +1; + } + feenox_call(feenox_problem_multifreedom_add(j_global, coefficients)); + feenox_free(coefficients); + } + } + +#endif + + return FEENOX_OK; +} diff --git a/src/pdes/neutron_sn/bulk.c b/src/pdes/neutron_sn/bulk.c new file mode 100644 index 00000000..d89527d2 --- /dev/null +++ b/src/pdes/neutron_sn/bulk.c @@ -0,0 +1,231 @@ +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: bulk elements + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#include "feenox.h" +#include "neutron_sn.h" + + +int feenox_problem_build_allocate_aux_neutron_sn(unsigned int n_nodes) { + + neutron_sn.n_nodes = n_nodes; + int dofs = neutron_sn.groups * neutron_sn.directions; + int size = neutron_sn.n_nodes * neutron_sn.groups * neutron_sn.directions; + + if (neutron_sn.Ki != NULL) { + gsl_matrix_free(neutron_sn.Ki); + gsl_matrix_free(neutron_sn.Ai); + gsl_matrix_free(neutron_sn.Xi); + if (neutron_sn.has_sources) { + gsl_vector_free(neutron_sn.Si); + } + + + gsl_matrix_free(neutron_sn.P); + gsl_matrix_free(neutron_sn.OMEGAB); + gsl_matrix_free(neutron_sn.AH); + gsl_matrix_free(neutron_sn.Xi); + if (neutron_sn.has_fission) { + gsl_matrix_free(neutron_sn.XH); + } + } + + feenox_check_alloc(neutron_sn.Ki = gsl_matrix_calloc(size, size)); + feenox_check_alloc(neutron_sn.Ai = gsl_matrix_calloc(size, size)); + feenox_check_alloc(neutron_sn.Xi = gsl_matrix_calloc(size, size)); + if (neutron_sn.has_sources) { + feenox_check_alloc(neutron_sn.Si = gsl_vector_calloc(size)); + } + + feenox_check_alloc(neutron_sn.P = gsl_matrix_calloc(dofs, size)); + feenox_check_alloc(neutron_sn.OMEGAB = gsl_matrix_calloc(dofs, size)); + feenox_check_alloc(neutron_sn.AH = gsl_matrix_calloc(dofs, size)); + if (neutron_sn.has_fission) { + feenox_check_alloc(neutron_sn.XH = gsl_matrix_calloc(dofs, size)); + } + + return FEENOX_OK; +} + +int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigned int q) { + +#ifdef HAVE_PETSC + + feenox_call(feenox_mesh_compute_wHB_at_gauss(e, q)); + double *x = feenox_mesh_compute_x_at_gauss_if_needed(e, q, neutron_sn.space_XS); + material_t *material = feenox_mesh_get_material(e); + + // initialize to zero + gsl_matrix_set_zero(neutron_sn.removal); + if (neutron_sn.has_fission) { + gsl_matrix_set_zero(neutron_sn.nufission); + } + if (neutron_sn.has_sources) { + gsl_vector_set_zero(neutron_sn.src); + } + + // TODO: if XS are uniform then compute only once these cached properties + for (int g = 0; g < neutron_sn.groups; g++) { + // independent sources + if (neutron_sn.S[g].defined) { + neutron_sn.S[g].eval(&neutron_sn.S[g], x, material); + } + + // fission + if (neutron_sn.nu_Sigma_f[g].defined) { + neutron_sn.nu_Sigma_f[g].eval(&neutron_sn.nu_Sigma_f[g], x, material); + } + + // scattering + for (int g_prime = 0; g_prime < neutron_sn.groups; g_prime++) { + if (neutron_sn.Sigma_s0[g_prime][g].defined) { + neutron_sn.Sigma_s0[g_prime][g].eval(&neutron_sn.Sigma_s0[g_prime][g], x, material); + } + if (neutron_sn.Sigma_s1[g_prime][g].defined) { + neutron_sn.Sigma_s1[g_prime][g].eval(&neutron_sn.Sigma_s1[g_prime][g], x, material); + } + } + + // absorption + if (neutron_sn.Sigma_t[g].defined) { + neutron_sn.Sigma_t[g].eval(&neutron_sn.Sigma_t[g], x, material); + } + if (neutron_sn.Sigma_a[g].defined) { + neutron_sn.Sigma_a[g].eval(&neutron_sn.Sigma_a[g], x, material); + } + } + + for (unsigned int m = 0; m < neutron_sn.directions; m++) { + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + // independent sources + int diag = dof_index(m,g); + gsl_vector_set(neutron_sn.src, diag, neutron_sn.S[g].value); + + // scattering and fission + for (unsigned int g_prime = 0; g_prime < neutron_sn.groups; g_prime++) { + for (unsigned int n_prime = 0; n_prime < neutron_sn.directions; n_prime++) { + // scattering + double xi = -neutron_sn.w[n_prime] * neutron_sn.Sigma_s0[g_prime][g].value; + // anisotropic scattering l = 1 + if (neutron_sn.Sigma_s1[g_prime][g].defined) { + xi -= neutron_sn.w[n_prime] * neutron_sn.Sigma_s1[g_prime][g].value * 3.0 * feenox_mesh_dot(neutron_sn.Omega[m], neutron_sn.Omega[n_prime]); + } + gsl_matrix_set(neutron_sn.removal, diag, dof_index(n_prime,g_prime), xi); + + // fision + if (neutron_sn.has_fission) { + gsl_matrix_set(neutron_sn.nufission, diag, dof_index(n_prime,g_prime), +neutron_sn.w[n_prime] * neutron_sn.chi[g] * neutron_sn.nu_Sigma_f[g_prime].value); + } + } + } + + // absorption + double xi = gsl_matrix_get(neutron_sn.removal, diag, diag); + if (neutron_sn.Sigma_t[g].defined) { + xi += neutron_sn.Sigma_t[g].value; + } else { + xi += neutron_sn.Sigma_a[g].value; + for (unsigned int g_prime = 0; g_prime < neutron_sn.groups; g_prime++) { + xi += neutron_sn.Sigma_s0[g][g_prime].value; + } + } + gsl_matrix_set(neutron_sn.removal, diag, diag, xi); + } + } + + // TODO: store current number of nodes in feenox.pde and call the method automatically? + // this should be done outside the gauss loop! + if (neutron_sn.n_nodes != e->type->nodes) { + feenox_call(feenox_problem_build_allocate_aux_neutron_sn(e->type->nodes)); + } + + // initialize Ki Ai Xi Si <- 0 + gsl_matrix_set_zero(neutron_sn.Ki); + gsl_matrix_set_zero(neutron_sn.Ai); + if (neutron_sn.has_fission) { + gsl_matrix_set_zero(neutron_sn.Xi); + } + if (neutron_sn.has_sources) { + gsl_vector_set_zero(neutron_sn.Si); + } + + // petrov stabilization + int MG = neutron_sn.directions * neutron_sn.groups; +// double tau = 0.5 * feenox_var_value(neutron_sn.sn_alpha) * neutron_sn.supg_dimension_factor * e->type->volume(e) / e->type->area(e); + double tau = feenox_var_value(neutron_sn.sn_alpha) * e->type->size(e); + + feenox_call(gsl_matrix_memcpy(neutron_sn.P, e->type->H_Gc[q])); + for (unsigned int j = 0; j < neutron_sn.n_nodes; j++) { + int NGj = MG*j; + for (unsigned int n = 0; n < neutron_sn.directions; n++) { + for (unsigned int m = 0; m < feenox.pde.dim; m++) { + double xi = tau * neutron_sn.Omega[n][m] * gsl_matrix_get(e->B[q], m, j); + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + int diag = dof_index(n,g); + gsl_matrix_add_to_element(neutron_sn.P, diag, NGj + diag, xi); + } + } + } + } + + // petrov-stabilized leakage term + feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_sn.OMEGA, e->B_G[q], 0.0, neutron_sn.OMEGAB)); + feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, e->w[q], neutron_sn.P, neutron_sn.OMEGAB, 1.0, neutron_sn.Ki)); + + // petrov-stabilized revmoval term + feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_sn.removal, e->type->H_Gc[q], 0.0, neutron_sn.AH)); + feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, e->w[q], neutron_sn.P, neutron_sn.AH, 1.0, neutron_sn.Ai)); + + // petrov-stabilized fission term + if (neutron_sn.has_fission) { + feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_sn.nufission, e->type->H_Gc[q], 0.0, neutron_sn.XH)); + feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, e->w[q], neutron_sn.P, neutron_sn.XH, 1.0, neutron_sn.Xi)); + } + // petrov-stabilized source term + if (neutron_sn.has_sources) { + feenox_call(gsl_blas_dgemv(CblasTrans, e->w[q], neutron_sn.P, neutron_sn.src, 1.0, neutron_sn.Si)); + } + + // for source-driven problems + // K = Ki + Ai - Xi + // for criticallity problems + // K = Ki + Ai + // M = Xi + gsl_matrix_add(neutron_sn.Ki, neutron_sn.Ai); + if (neutron_sn.has_fission) { + if (neutron_sn.has_sources) { + gsl_matrix_scale(neutron_sn.Xi, -1.0); + gsl_matrix_add(neutron_sn.Ki, neutron_sn.Xi); + } else { + gsl_matrix_add(feenox.pde.Mi, neutron_sn.Xi); + } + } + + if (neutron_sn.has_sources) { + gsl_vector_add(feenox.pde.bi, neutron_sn.Si); + } + + gsl_matrix_add(feenox.pde.Ki, neutron_sn.Ki); +#endif + + return FEENOX_OK; + +} + diff --git a/src/pdes/neutron_sn/init.c b/src/pdes/neutron_sn/init.c new file mode 100644 index 00000000..b23b081a --- /dev/null +++ b/src/pdes/neutron_sn/init.c @@ -0,0 +1,566 @@ +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: initialization + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#include "feenox.h" +#include "neutron_sn.h" +neutron_sn_t neutron_sn; + +int feenox_problem_neutron_sn_init_cosines(double mu1) { + // equation 21 in stammler + // equation 4-8 in lewis + int N_over_2 = neutron_sn.N/2; + double *mu = NULL; + feenox_check_alloc(mu = calloc(N_over_2, sizeof(double))); + + mu[0] = mu1; + double C = 2*(1-3*gsl_pow_2(mu1))/(neutron_sn.N-2); + for (int i = 1; i < N_over_2; i++) { + mu[i] = sqrt(gsl_pow_2(mu1) + C*i); + } + + // we have these triangles + // S4 + // 1 + // 1 1 + // + // S6 + // 1 + // 2 2 + // 1 2 1 + // + // S8 + // 1 + // 2 2 + // 2 3 2 + // 1 2 2 1 + + int i = 0; + int j = 0; + int k = 0; + int n = 0; + for (int row = 0; row < N_over_2; row++) { + for (int col = 0; col <= row; col++) { + i = N_over_2 - row - 1; + j = col; + k = N_over_2 + 2 - 3 - i - j; + +// printf("%d\t%d %d %d\n", n, i, j, k); + + neutron_sn.Omega[n][0] = mu[i]; + neutron_sn.Omega[n][1] = mu[j]; + neutron_sn.Omega[n][2] = (feenox.pde.dim == 3) ? mu[k] : 0; + n++; + + } + } + + feenox_free(mu); + + return FEENOX_OK; +} + +int feenox_problem_init_parser_neutron_sn(void) { + +///kw_pde+PROBLEM+detail * `neutron_sn` multi-group core-level neutron transport using +///kw_pde+PROBLEM+detail - discrete ordinates $S_N$ for angular discretization, and +///kw_pde+PROBLEM+detail - isoparametric finite elements for spatial discretization. + + // we are FEM + feenox.mesh.default_field_location = field_location_nodes; + +#ifdef HAVE_PETSC + feenox.pde.init_runtime_particular = feenox_problem_init_runtime_neutron_sn; + feenox.pde.bc_parse = feenox_problem_bc_parse_neutron_sn; +#ifdef HAVE_SLEPC + feenox.pde.setup_eps = feenox_problem_setup_eps_neutron_sn; +#endif + feenox.pde.setup_ksp = feenox_problem_setup_ksp_neutron_sn; + feenox.pde.setup_pc = feenox_problem_setup_pc_neutron_sn; + feenox.pde.build_element_volumetric_gauss_point = feenox_problem_build_volumetric_gauss_point_neutron_sn; + feenox.pde.solve_post = feenox_problem_solve_post_neutron_sn; + + // default is 1 group + if (neutron_sn.groups == 0) { + neutron_sn.groups = 1; + } + + // default is N=2 + if (neutron_sn.N == 0) { + neutron_sn.N = 2; + } + if (neutron_sn.N % 2 != 0) { + feenox_push_error_message("number of ordinates N = %d has to be even", neutron_sn.N); + return FEENOX_ERROR; + } + + // stammler's eq 6.19 + switch(feenox.pde.dim) { + case 1: + neutron_sn.directions = neutron_sn.N; +// neutron_sn.supg_dimension_factor = 1; + break; + case 2: + neutron_sn.directions = 0.5*neutron_sn.N*(neutron_sn.N+2); +// neutron_sn.supg_dimension_factor = 4; +// neutron_sn.supg_dimension_factor = 6; + break; + case 3: + neutron_sn.directions = neutron_sn.N*(neutron_sn.N+2); +// neutron_sn.supg_dimension_factor = 6; +// neutron_sn.supg_dimension_factor = 12; +// neutron_sn.supg_dimension_factor = 16; + break; + } + + + // dofs = number of directions * number of groups + feenox.pde.dofs = neutron_sn.directions * neutron_sn.groups; + + // ------- neutron transport outputs ----------------------------------- + + // the angular fluxes psi + feenox_check_alloc(feenox.pde.unknown_name = calloc(feenox.pde.dofs, sizeof(char *))); + // TODO: document from comments + for (int n = 0; n < neutron_sn.directions; n++) { + for (int g = 0; g < neutron_sn.groups; g++) { + feenox_check_minusone(asprintf(&feenox.pde.unknown_name[n * neutron_sn.groups + g], "psi%d.%d", n+1, g+1)); + } + } + feenox_call(feenox_problem_define_solutions()); + + // the scalar fluxes psi + // TODO: document + // TODO: for one group make an alias between phi1 and phi + feenox_check_alloc(neutron_sn.phi = calloc(neutron_sn.groups, sizeof(function_t *))); + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + char *name = NULL; + feenox_check_minusone(asprintf(&name, "phi%u", g+1)); + feenox_call(feenox_problem_define_solution_function(name, &neutron_sn.phi[g], FEENOX_SOLUTION_NOT_GRADIENT)); + feenox_free(name); + } + +///va_neutron_sn+keff+desc The effective multiplication factor\ $k_\text{eff}$. + neutron_sn.keff = feenox_define_variable_get_ptr("keff"); + +///va_neutron_sn+sn_alpha+desc The stabilization parameter\ $\alpha$ for $S_N$. + neutron_sn.sn_alpha = feenox_define_variable_get_ptr("sn_alpha"); + feenox_var_value(neutron_sn.sn_alpha) = 0.5; + + // ------------ initialize the SN weights ------------------------------------ + // TODO: when adding FVM, these are the same so we should move it into a single location + feenox_check_alloc(neutron_sn.w = calloc(neutron_sn.directions, sizeof(double))); + feenox_check_alloc(neutron_sn.Omega = calloc(neutron_sn.directions, sizeof(double *))); + for (int m = 0; m < neutron_sn.directions; m++) { + feenox_check_alloc(neutron_sn.Omega[m] = calloc(3, sizeof(double))); + } + + if (feenox.pde.dim == 1) { + // table 3-1 page 121 lewiss + // in 1D the directions are the zeros of the legendre polymonials + // and the weights are such the gauss quadrature is exact + switch (neutron_sn.N) { + case 2: + neutron_sn.Omega[0][0] = +1.0/M_SQRT3; + neutron_sn.w[0] = 1.0/2.0; + + neutron_sn.Omega[1][0] = -neutron_sn.Omega[0][0]; + neutron_sn.w[1] = neutron_sn.w[0]; + break; + case 4: + neutron_sn.Omega[0][0] = sqrt(3.0/7.0 - 2.0/7.0 * sqrt(6.0/5.0)); + neutron_sn.w[0] = 0.6521451549/2.0; + + neutron_sn.Omega[1][0] = sqrt(3.0/7.0 + 2.0/7.0 * sqrt(6.0/5.0)); + neutron_sn.w[1] = 0.3478548451/2.0; + + neutron_sn.Omega[2][0] = -neutron_sn.Omega[0][0]; + neutron_sn.w[2] = neutron_sn.w[0]; + + neutron_sn.Omega[3][0] = -neutron_sn.Omega[1][0]; + neutron_sn.w[3] = neutron_sn.w[1]; + break; + + case 6: + neutron_sn.Omega[0][0] = 0.2386191860; + neutron_sn.w[0] = 0.4679139346/2.0; + + neutron_sn.Omega[1][0] = 0.6612093864; + neutron_sn.w[1] = 0.3607615730/2.0; + + neutron_sn.Omega[2][0] = 0.9324695142; + neutron_sn.w[2] = 0.1713244924/2.0; + + neutron_sn.Omega[3][0] = -neutron_sn.Omega[0][0]; + neutron_sn.w[3] = neutron_sn.w[0]; + + neutron_sn.Omega[4][0] = -neutron_sn.Omega[1][0]; + neutron_sn.w[4] = neutron_sn.w[1]; + + neutron_sn.Omega[5][0] = -neutron_sn.Omega[2][0]; + neutron_sn.w[5] = neutron_sn.w[2]; + break; + + case 8: + neutron_sn.Omega[0][0] = 0.1834346424; + neutron_sn.w[0] = 0.3626837834/2.0; + + neutron_sn.Omega[1][0] = 0.5255324099; + neutron_sn.w[1] = 0.3137066459/2.0; + + neutron_sn.Omega[2][0] = 0.7966664774; + neutron_sn.w[2] = 0.2223810344/2.0; + + neutron_sn.Omega[3][0] = 0.9602898564; + neutron_sn.w[3] = 0.1012285363/2.0; + + neutron_sn.Omega[4][0] = -neutron_sn.Omega[0][0]; + neutron_sn.w[4] = neutron_sn.w[0]; + + neutron_sn.Omega[5][0] = -neutron_sn.Omega[1][0]; + neutron_sn.w[5] = neutron_sn.w[1]; + + neutron_sn.Omega[6][0] = -neutron_sn.Omega[2][0]; + neutron_sn.w[6] = neutron_sn.w[2]; + + neutron_sn.Omega[7][0] = -neutron_sn.Omega[3][0]; + neutron_sn.w[7] = neutron_sn.w[3]; + break; + + default: + feenox_push_error_message("unsupported N = %d", neutron_sn.N); + return FEENOX_ERROR; + break; + } + } else { + // table 4-1 page 162 from lewiss + // we could have used file sndir2.dat from fentraco, but the values are different + + // this is for 2 and 3 dimensions + int N_octs = (feenox.pde.dim == 2) ? 4 : 8; + int J_octs = neutron_sn.directions / N_octs; + + // first set the weights as all the possible permutations in the first octant + // and then we fill in the others + // the initial cosine directions passed to feenox_problem_neutron_sn_init_cosines() + // are taken from table 4-1 from lewiss (p 162) + // which coincide with table 1 p 208 from stammler-abbate + + switch (neutron_sn.N) { + case 2: + { + feenox_call(feenox_problem_neutron_sn_init_cosines(1.0/M_SQRT3)); + + neutron_sn.w[0] = 1.0/(double)(N_octs); + } + break; + + case 4: + { + feenox_call(feenox_problem_neutron_sn_init_cosines(0.3500212)); + double s4_w1 = 1.0/3.0; + + neutron_sn.w[0] = s4_w1/(double)(N_octs); + + neutron_sn.w[1] = s4_w1/(double)(N_octs); + neutron_sn.w[2] = s4_w1/(double)(N_octs); + } + break; + + case 6: + { + feenox_call(feenox_problem_neutron_sn_init_cosines(0.2666355)); + + double s6_w1 = 0.1761263; + double s6_w2 = 0.1572071; + + neutron_sn.w[0] = s6_w1/(double)(N_octs); + + neutron_sn.w[1] = s6_w2/(double)(N_octs); + neutron_sn.w[2] = s6_w2/(double)(N_octs); + + neutron_sn.w[3] = s6_w1/(double)(N_octs); + neutron_sn.w[4] = s6_w2/(double)(N_octs); + neutron_sn.w[5] = s6_w1/(double)(N_octs); + } + break; + + case 8: + { + feenox_call(feenox_problem_neutron_sn_init_cosines(0.2182179)); + + double s8_w1 = 0.1209877; + double s8_w2 = 0.0907407; + double s8_w3 = 0.0925926; + + + neutron_sn.w[0] = s8_w1/(double)(N_octs); + + neutron_sn.w[1] = s8_w2/(double)(N_octs); + neutron_sn.w[2] = s8_w2/(double)(N_octs); + + neutron_sn.w[3] = s8_w2/(double)(N_octs); + neutron_sn.w[4] = s8_w3/(double)(N_octs); + neutron_sn.w[5] = s8_w2/(double)(N_octs); + + neutron_sn.w[6] = s8_w1/(double)(N_octs); + neutron_sn.w[7] = s8_w2/(double)(N_octs); + neutron_sn.w[8] = s8_w2/(double)(N_octs); + neutron_sn.w[9] = s8_w1/(double)(N_octs); + + } + break; + + default: + feenox_push_error_message("unsupported N = %d in %dD", neutron_sn.N, feenox.pde.dim); + return FEENOX_ERROR; + break; + } + + + // we have filled the first octant, now fill in the other ones + for (int n = 1; n < N_octs; n++) { + for (int j = 0; j < J_octs; j++) { + neutron_sn.Omega[n*J_octs + j][0] = ((n & 1) ? (-1) : (+1)) * neutron_sn.Omega[j][0]; + neutron_sn.Omega[n*J_octs + j][1] = ((n & 2) ? (-1) : (+1)) * neutron_sn.Omega[j][1]; + neutron_sn.Omega[n*J_octs + j][2] = ((n & 4) ? (-1) : (+1)) * neutron_sn.Omega[j][2]; + neutron_sn.w[n*J_octs + j] = neutron_sn.w[j]; + } + } + } + + // checks + double s = 0; + for (unsigned int n = 0; n < neutron_sn.directions; n++) { + s += neutron_sn.w[n]; + } + assert(fabs(s-1) < 1e-6); + + for (unsigned int d = 0; d < feenox.pde.dim; d++) { + s = 0; + for (unsigned int n = 0; n < neutron_sn.directions; n++) { + s += neutron_sn.w[n] * neutron_sn.Omega[n][d]; + } + assert(fabs(s) < 1e-6); + + s = 0; + for (unsigned int n = 0; n < neutron_sn.directions; n++) { + s += neutron_sn.w[n] * gsl_pow_2(neutron_sn.Omega[n][d]); + } + assert(fabs(s-1.0/3.0) < 1e-6); + } + + // --------------------------------------------------------------------------- + +#endif + + return FEENOX_OK; +} + + +int feenox_problem_init_runtime_neutron_sn(void) { + +#ifdef HAVE_PETSC + feenox.pde.mesh->data_type = data_type_node; + feenox.pde.spatial_unknowns = feenox.pde.mesh->n_nodes; + feenox.pde.size_global = feenox.pde.spatial_unknowns * feenox.pde.dofs; + + // set the size of the eigenvectors (we did not know their size in init_parser() above + for (PetscInt i = 0; i < feenox.pde.nev; i++) { + feenox_call(feenox_vector_set_size(feenox.pde.vectors.phi[i], feenox.pde.size_global)); + } + + // initialize XSs + // TODO: allow not to give the group number in one-group problems + unsigned int G = neutron_sn.groups; + feenox_check_alloc(neutron_sn.Sigma_t = calloc(G, sizeof(distribution_t))); + feenox_check_alloc(neutron_sn.Sigma_a = calloc(G, sizeof(distribution_t))); + feenox_check_alloc(neutron_sn.nu_Sigma_f = calloc(G, sizeof(distribution_t))); + feenox_check_alloc(neutron_sn.S = calloc(G, sizeof(distribution_t))); + feenox_check_alloc(neutron_sn.Sigma_s0 = calloc(G, sizeof(distribution_t *))); + feenox_check_alloc(neutron_sn.Sigma_s1 = calloc(G, sizeof(distribution_t *))); + for (unsigned int g = 0; g < G; g++) { + char *name = NULL; + + feenox_check_minusone(asprintf(&name, "Sigma_t%d", g+1)); + feenox_distribution_init(&neutron_sn.Sigma_t[g], name); + if (neutron_sn.Sigma_t[g].defined) { + neutron_sn.Sigma_t[g].uniform = feenox_expression_depends_on_space(neutron_sn.Sigma_t[g].dependency_variables); + } + feenox_free(name); + + feenox_check_minusone(asprintf(&name, "Sigma_a%d", g+1)); + feenox_distribution_init(&neutron_sn.Sigma_a[g], name); + if (neutron_sn.Sigma_a[g].defined) { + neutron_sn.Sigma_a[g].uniform = feenox_expression_depends_on_space(neutron_sn.Sigma_a[g].dependency_variables); + } + feenox_free(name); + + feenox_check_minusone(asprintf(&name, "nuSigma_f%d", g+1)); + feenox_distribution_init(&neutron_sn.nu_Sigma_f[g], name); + if (neutron_sn.nu_Sigma_f[g].defined) { + neutron_sn.has_fission = 1; + neutron_sn.nu_Sigma_f[g].uniform = feenox_expression_depends_on_space(neutron_sn.nu_Sigma_f[g].dependency_variables); + } + feenox_free(name); + + feenox_check_minusone(asprintf(&name, "S%d", g+1)); + feenox_distribution_init(&neutron_sn.S[g], name); + if (neutron_sn.S[g].defined) { + neutron_sn.has_sources = 1; + neutron_sn.S[g].uniform = feenox_expression_depends_on_space(neutron_sn.S[g].dependency_variables); + } + feenox_free(name); + + feenox_check_alloc(neutron_sn.Sigma_s0[g] = calloc(G, sizeof(distribution_t))); + feenox_check_alloc(neutron_sn.Sigma_s1[g] = calloc(G, sizeof(distribution_t))); + for (int g_prime = 0; g_prime < G; g_prime++) { + feenox_check_minusone(asprintf(&name, "Sigma_s%d.%d", g+1, g_prime+1)); + feenox_distribution_init(&neutron_sn.Sigma_s0[g][g_prime], name); + if (neutron_sn.Sigma_s0[g][g_prime].defined) { + neutron_sn.Sigma_s0[g][g_prime].uniform = feenox_expression_depends_on_space(neutron_sn.Sigma_s0[g][g_prime].dependency_variables); + } + feenox_free(name); + + feenox_check_minusone(asprintf(&name, "Sigma_s_one%d.%d", g+1, g_prime+1)); + feenox_distribution_init(&neutron_sn.Sigma_s1[g][g_prime], name); + if (neutron_sn.Sigma_s1[g][g_prime].defined) { + neutron_sn.Sigma_s1[g][g_prime].uniform = feenox_expression_depends_on_space(neutron_sn.Sigma_s1[g][g_prime].dependency_variables); + } + feenox_free(name); + + } + } + + + if (neutron_sn.has_sources == 0 && neutron_sn.has_fission == 0) { + feenox_push_error_message("neither fission nor sources found"); + return FEENOX_ERROR; + } + + if (neutron_sn.has_sources == 0) { +#ifdef HAVE_SLEPC + feenox.pde.math_type = math_type_eigen; + feenox.pde.solve = feenox_problem_solve_slepc_eigen; + feenox.pde.has_mass = 1; + feenox.pde.has_rhs = 0; + + // TODO: higher harmonics? + feenox.pde.nev = 1; + + // define eigenvectors (we don't know its size yet) + feenox_check_alloc(feenox.pde.vectors.phi = calloc(feenox.pde.nev, sizeof(vector_t *))); + for (int g = 0; g < feenox.pde.nev; g++) { + char *modename = NULL; + feenox_check_minusone(asprintf(&modename, "eig%d", g+1)); + feenox_check_alloc(feenox.pde.vectors.phi[g] = feenox_define_vector_get_ptr(modename, 0)); + feenox_free(modename); + } +#else + feenox_push_error_message("criticality problems cannot be solved without SLEPc"); + return FEENOX_ERROR; +#endif + } else { + // TODO: non-linear? + feenox.pde.math_type = math_type_linear; + feenox.pde.solve = feenox_problem_solve_petsc_linear; + feenox.pde.has_mass = 0; + feenox.pde.has_rhs = 1; + } + + // allocate elemental XS matrices + unsigned int N = neutron_sn.directions; + unsigned int NG = N*G; + feenox_check_alloc(neutron_sn.removal = gsl_matrix_calloc(NG, NG)); + feenox_check_alloc(neutron_sn.nufission = gsl_matrix_calloc(NG, NG)); + feenox_check_alloc(neutron_sn.src = gsl_vector_calloc(NG)); + feenox_check_alloc(neutron_sn.chi = calloc(G, sizeof(double))); + // TODO: read a vector called "chi" + neutron_sn.chi[0] = 1; + + // direction matrix + feenox_check_alloc(neutron_sn.OMEGA = gsl_matrix_calloc(NG, NG * feenox.pde.dim)); + for (unsigned int n = 0; n < N; n++) { + for (unsigned int g = 0; g < G; g++) { + for (unsigned int m = 0; m < feenox.pde.dim; m++) { + gsl_matrix_set(neutron_sn.OMEGA, dof_index(n,g), m*NG + dof_index(n,g), neutron_sn.Omega[n][m]); + } + } + } + + + + feenox.pde.has_stiffness = 1; + + feenox.pde.has_jacobian_K = 0; + feenox.pde.has_jacobian_M = 0; + feenox.pde.has_jacobian_b = 0; + feenox.pde.has_jacobian = feenox.pde.has_jacobian_K || feenox.pde.has_jacobian_M || feenox.pde.has_jacobian_b; + + feenox.pde.symmetric_K = 0; + feenox.pde.symmetric_M = 0; +#endif + return FEENOX_OK; +} + +#ifdef HAVE_PETSC +int feenox_problem_setup_pc_neutron_sn(PC pc) { + + PCType pc_type = NULL; + petsc_call(PCGetType(pc, &pc_type)); + if (pc_type == NULL) { + petsc_call(PCSetType(pc, PCLU)); +#ifdef PETSC_HAVE_MUMPS + petsc_call(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); +#endif + } + + return FEENOX_OK; +} + +int feenox_problem_setup_ksp_neutron_sn(KSP ksp ) { + + KSPType ksp_type = NULL; + petsc_call(KSPGetType(ksp, &ksp_type)); + if (ksp_type == NULL) { + petsc_call(KSPSetType(ksp, KSPPREONLY)) + } + + return FEENOX_OK; +} +#endif + +#ifdef HAVE_SLEPC +int feenox_problem_setup_eps_neutron_sn(EPS eps) { + + // generalized non-hermitian problem +// petsc_call(EPSSetProblemType(eps, EPS_GNHEP)); + + if (feenox.pde.eigen_formulation == eigen_formulation_omega) { + petsc_call(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE)); + } else { + petsc_call(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE)); + } + + + return FEENOX_OK; +} +#endif + diff --git a/src/pdes/neutron_sn/methods.h b/src/pdes/neutron_sn/methods.h new file mode 100644 index 00000000..10f62529 --- /dev/null +++ b/src/pdes/neutron_sn/methods.h @@ -0,0 +1,52 @@ +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: virtual methods + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#ifndef NEUTRON_SN_METHODS_H +#define NEUTRON_SN_METHODS_H + +// neutron_sn/parser.c +extern int feenox_problem_parse_problem_neutron_sn(const char *); +extern int feenox_problem_parse_post_neutron_sn(mesh_write_t *mesh_write, const char *token); + +// neutron_sn/init.c +extern int feenox_problem_neutron_sn_init_cosines(double mu1); +extern int feenox_problem_init_parser_neutron_sn(void); +extern int feenox_problem_init_runtime_neutron_sn(void); +#ifdef HAVE_PETSC +extern int feenox_problem_setup_pc_neutron_sn(PC pc); +extern int feenox_problem_setup_ksp_neutron_sn(KSP ksp); +#endif +#ifdef HAVE_SLEPC +extern int feenox_problem_setup_eps_neutron_sn(EPS eps); +#endif + +// neutron_sn/bulk.c +extern int feenox_problem_build_allocate_aux_neutron_sn(unsigned int n_nodes); +extern int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *element, unsigned int q); + +// neutron_sn/bc.c +extern int feenox_problem_bc_parse_neutron_sn(bc_data_t *bc_data, const char *lhs, char *rhs); +extern int feenox_problem_bc_set_neutron_sn_vacuum(bc_data_t *bc_data, element_t *e, size_t j_global); +extern int feenox_problem_bc_set_neutron_sn_mirror(bc_data_t *bc_data, element_t *e, size_t j_global); + +// neutron_sn/post.c +extern int feenox_problem_solve_post_neutron_sn(void); +#endif // NEUTRON_SN_METHODS_H diff --git a/src/pdes/neutron_sn/neutron_sn.h b/src/pdes/neutron_sn/neutron_sn.h new file mode 100644 index 00000000..b749751f --- /dev/null +++ b/src/pdes/neutron_sn/neutron_sn.h @@ -0,0 +1,115 @@ +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: global header + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#ifndef NEUTRON_SN_H +#define NEUTRON_SN_H + +#define dof_index(n,g) ((n)*neutron_sn.groups + (g)) + +typedef struct neutron_sn_t neutron_sn_t; + + + +struct neutron_sn_t { + + unsigned int N; + unsigned int directions; + unsigned int groups; + + // directions and weights + double **Omega; + double *w; + + // for SUPG: factor relating volume/area of an element depending on the dimension + // if there were circles, spheres: 1 for 1D, 4 for 2D and 6 for 3D + // if there were triang, tets: 1 for 1D, 6 for 2D and 12 for 3d +// double supg_dimension_factor; + + + // total XS + distribution_t *Sigma_t; + + // absorption XS + distribution_t *Sigma_a; + + // scattering XS + distribution_t **Sigma_s0; + distribution_t **Sigma_s1; + + // product of nu and fission XS + distribution_t *nu_Sigma_f; + + // product of energy per fission and fission XS +// distribution_t *e_Sigma_f; + + // independent neutron source + distribution_t *S; + + + // removal XSs: groups x groups + gsl_matrix *removal; + + // fission XSs: groups x groups + gsl_matrix *nufission; + + // independent sources: groups + gsl_vector *src; + + // fission spectrum: groups + double *chi; + + + + // elemental matrices + unsigned int n_nodes; + + gsl_matrix *P; + gsl_matrix *OMEGA; + gsl_matrix *OMEGAB; + gsl_matrix *AH; + gsl_matrix *XH; + + gsl_matrix *Ki; + gsl_matrix *Ai; + gsl_matrix *Xi; + gsl_vector *Si; + + + int has_sources; + int has_fission; + int space_XS; + + var_t *sn_alpha; + + // --- results ---------------------- + + // effective multiplication factor + var_t *keff; + // scalar fluxes (array of G functions) + function_t **phi; + + +}; +extern neutron_sn_t neutron_sn; + + +#endif /* NEUTRON_SN_H */ + diff --git a/src/pdes/neutron_sn/parser.c b/src/pdes/neutron_sn/parser.c new file mode 100644 index 00000000..4749963d --- /dev/null +++ b/src/pdes/neutron_sn/parser.c @@ -0,0 +1,85 @@ +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: input parsing + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#include "feenox.h" +#include "neutron_sn.h" +#include "../../parser/parser.h" + +int feenox_problem_parse_problem_neutron_sn(const char *token) { + + if (token != NULL) { + if (strcasecmp(token, "GROUPS") == 0) { + double xi = 0; + feenox_call(feenox_parser_expression_in_string(&xi)); + neutron_sn.groups = (unsigned int)(xi); + + } else if (strcasecmp(token, "SN") == 0 || strcasecmp(token, "N") == 0) { + double xi = 0; + feenox_call(feenox_parser_expression_in_string(&xi)); + neutron_sn.N = (unsigned int)(xi); + if (neutron_sn.N != 2 && + neutron_sn.N != 4 && + neutron_sn.N != 6 && + neutron_sn.N != 8) { + feenox_push_error_message("S%d is not supported, only S2, S4, S6 and S8 are", neutron_sn.N); + return FEENOX_ERROR; + } + + } else { + feenox_push_error_message("undefined keyword '%s'", token); + return FEENOX_ERROR; + } + } + + return FEENOX_OK; +} + + + +int feenox_problem_parse_post_neutron_sn(mesh_write_t *mesh_write, const char *token) { + + if (strcmp(token, "all") == 0) { + feenox_call(feenox_problem_parse_post_neutron_sn(mesh_write, "scalar_fluxes")); + feenox_call(feenox_problem_parse_post_neutron_sn(mesh_write, "angular_fluxes")); + + } else if (strcmp(token, "scalar_flux") == 0 || strcmp(token, "scalar_fluxes") == 0) { + + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + feenox_call(feenox_add_post_field(mesh_write, 1, &neutron_sn.phi[g]->name, NULL, field_location_nodes)); + } + + } else if (strcmp(token, "angular_flux") == 0 || strcmp(token, "angular_fluxes") == 0) { + + for (int n = 0; n < neutron_sn.directions; n++) { + for (int g = 0; g < neutron_sn.groups; g++) { + feenox_call(feenox_add_post_field(mesh_write, 1, &feenox.pde.unknown_name[n * neutron_sn.groups + g], NULL, field_location_nodes)); + } + } + + + } else { + feenox_push_error_message("undefined keyword '%s' for neutron_sn WRITE_RESULTS", token); + return FEENOX_ERROR; + } + + return FEENOX_OK; +} + diff --git a/src/pdes/neutron_sn/post.c b/src/pdes/neutron_sn/post.c new file mode 100644 index 00000000..da163356 --- /dev/null +++ b/src/pdes/neutron_sn/post.c @@ -0,0 +1,86 @@ + +/*------------ -------------- -------- --- ----- --- -- - - + * feenox's routines for neutron transport FEM: post + * + * Copyright (C) 2023 jeremy theler + * + * This file is part of FeenoX . + * + * feenox is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * FeenoX is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with FeenoX. If not, see . + *------------------- ------------ ---- -------- -- - - - + */ +#include "feenox.h" +#include "neutron_sn.h" + +int feenox_problem_solve_post_neutron_sn(void) { + +#ifdef HAVE_PETSC + if (neutron_sn.has_sources == 0) { +#ifdef HAVE_SLEPC + if (feenox.pde.eigen_formulation == eigen_formulation_omega) { + feenox_var_value(neutron_sn.keff) = 1.0/feenox.pde.eigenvalue[0]; + } else { + feenox_var_value(neutron_sn.keff) = feenox.pde.eigenvalue[0]; + } + + // TODO: this is wrong + // normalization without power + double num = 0; + double den = 0; + unsigned int g = 0; + size_t i = 0; + for (i = 0; i < feenox.pde.mesh->n_elements; i++) { + element_t *element = &feenox.pde.mesh->element[i]; + if (element->type != NULL && element->type->dim == feenox.pde.dim) { + num += element->type->volume(element); + for (g = 0; g < feenox.pde.dofs; g++) { + den += feenox_mesh_integral_over_element(element, feenox.pde.mesh, feenox.pde.solution[g]); + } + } + } + + double factor = num/den; + + // normalize the fluxes + for (size_t j = 0; j < feenox.pde.mesh->n_nodes; j++) { + for (g = 0; g < feenox.pde.dofs; g++) { + double xi = feenox_vector_get(feenox.pde.solution[g]->vector_value, j); + feenox_vector_set(feenox.pde.solution[g]->vector_value, j, factor*xi); + } + } +#endif + } + + + // compute the scalar fluxes out of the directional fluxes + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + feenox_problem_fill_aux_solution(neutron_sn.phi[g]); + } + + for (size_t j = 0; j < feenox.pde.mesh->n_nodes; j++) { + for (unsigned int g = 0; g < neutron_sn.groups; g++) { + double xi = 0; + for (unsigned int n = 0; n < neutron_sn.directions; n++) { + xi += neutron_sn.w[n] * feenox_vector_get(feenox.pde.solution[n * neutron_sn.groups + g]->vector_value, j); + } + // TODO: wrapper, the pde has to set the function value not the vector + feenox_vector_set(neutron_sn.phi[g]->vector_value, j, xi); + } + } + + +#endif + + return FEENOX_OK; +}