Skip to content

Commit

Permalink
renamed variables in sn
Browse files Browse the repository at this point in the history
  • Loading branch information
@gtheler
gtheler committed Jul 24, 2023
1 parent 5edd4e5 commit 275f76c
Show file tree
Hide file tree
Showing 7 changed files with 67 additions and 87 deletions.
30 changes: 0 additions & 30 deletions src/mesh/fem.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -607,36 +607,6 @@ int feenox_mesh_compute_B_at_gauss(element_t *e, unsigned int q, int integration
return FEENOX_OK;
}

/*
int feenox_mesh_compute_B_Gc_at_gauss(element_type_t *element_type, unsigned int q, int integration) {
unsigned int G = feenox.pde.dofs;
if (element_type->B_Gc == NULL) {
if (G == 1) {
element_type->B_Gc = element_type->B_c;
return FEENOX_OK;
}
feenox_check_alloc(element_type->B_Gc = calloc(element_type->gauss[integration].Q, sizeof(gsl_matrix *)));
}
unsigned int J = element_type->nodes;
unsigned int D = element_type->dim;
feenox_check_alloc(element_type->B_Gc[q] = gsl_matrix_calloc(G*D, G*D*J));
for (unsigned int d = 0; d < element_type->dim; d++) {
size_t Gd = G*d;
for (unsigned int j = 0; j < element_type->nodes; j++) {
size_t Gj = G*j;
double dhdxi = gsl_matrix_get(element_type->gauss[integration].B_c[q], d, j);
for (unsigned int g = 0; g < feenox.pde.dofs; g++) {
gsl_matrix_set(element_type->B_Gc[q], Gd+g, Gj+g, dhdxi);
}
}
}
return FEENOX_OK;
}
*/

int feenox_mesh_compute_B_G_at_gauss(element_t *e, unsigned int q, int integration) {

Expand Down
6 changes: 3 additions & 3 deletions src/pdes/neutron_diffusion/bulk.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_diffusion(element_t *e,
double *x = feenox_mesh_compute_x_at_gauss_if_needed(e, q, neutron_diffusion.space_XS);
material_t *material = feenox_mesh_get_material(e);

gsl_matrix_set_zero(neutron_diffusion.D_prime);
gsl_matrix_set_zero(neutron_diffusion.D_G);
gsl_matrix_set_zero(neutron_diffusion.R);
if (neutron_diffusion.has_fission) {
gsl_matrix_set_zero(neutron_diffusion.X);
Expand Down Expand Up @@ -108,7 +108,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_diffusion(element_t *e,
}

unsigned int index = m*neutron_diffusion.groups + g;
gsl_matrix_set(neutron_diffusion.D_prime, index, index, xi);
gsl_matrix_set(neutron_diffusion.D_G, index, index, xi);
}
}

Expand All @@ -124,7 +124,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_diffusion(element_t *e,

// elemental stiffness for the diffusion term B'*D*B
// TODO: convenience call for A'*B*A? that'd need an intermediate alloc
feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_diffusion.D_prime, e->B_G[q], 0.0, neutron_diffusion.DB));
feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_diffusion.D_G, e->B_G[q], 0.0, neutron_diffusion.DB));
feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, e->w[q], e->B_G[q], neutron_diffusion.DB, 1.0, neutron_diffusion.Li));

// elemental scattering H'*A*H
Expand Down
2 changes: 1 addition & 1 deletion src/pdes/neutron_diffusion/init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,7 @@ int feenox_problem_init_runtime_neutron_diffusion(void) {
}

// allocate elemental XS matrices
feenox_check_alloc(neutron_diffusion.D_prime = gsl_matrix_calloc(G * feenox.pde.dim, G * feenox.pde.dim));
feenox_check_alloc(neutron_diffusion.D_G = gsl_matrix_calloc(G * feenox.pde.dim, G * feenox.pde.dim));
feenox_check_alloc(neutron_diffusion.R = gsl_matrix_calloc(G, G));
feenox_check_alloc(neutron_diffusion.X = gsl_matrix_calloc(G, G));
feenox_check_alloc(neutron_diffusion.s = gsl_vector_calloc(G));
Expand Down
4 changes: 2 additions & 2 deletions src/pdes/neutron_diffusion/neutron_diffusion.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ struct neutron_diffusion_t {
distribution_t *S;

// diffusion coefficients, (groups * dim) x (groups * dim)
gsl_matrix *D_prime;
gsl_matrix *D_G;

// removal XSs: groups x groups
gsl_matrix *R;
Expand All @@ -51,7 +51,7 @@ struct neutron_diffusion_t {
double *chi;

unsigned int n_nodes;
// elemental matrices (size n_nodes*G x n_nodes*G)
// elemental matrices (size J*G x J*G)
gsl_matrix *Li; // leakage
gsl_matrix *Ai; // absorption
gsl_matrix *Fi; // fission
Expand Down
69 changes: 36 additions & 33 deletions src/pdes/neutron_sn/bulk.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,33 +32,36 @@ int feenox_problem_build_allocate_aux_neutron_sn(unsigned int J) {
int MG = M*G;
int JMG = J * MG;

if (neutron_sn.Ki != NULL) {
gsl_matrix_free(neutron_sn.Ki);
if (neutron_sn.Li != NULL) {
gsl_matrix_free(neutron_sn.Li);
gsl_matrix_free(neutron_sn.Ai);
gsl_matrix_free(neutron_sn.Xi);
if (neutron_sn.has_fission) {
gsl_matrix_free(neutron_sn.Fi);
}
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.DB);
gsl_matrix_free(neutron_sn.AH);
gsl_matrix_free(neutron_sn.Xi);
if (neutron_sn.has_fission) {
gsl_matrix_free(neutron_sn.Fi);
gsl_matrix_free(neutron_sn.XH);
}
}

feenox_check_alloc(neutron_sn.Ki = gsl_matrix_calloc(JMG, JMG));
feenox_check_alloc(neutron_sn.Li = gsl_matrix_calloc(JMG, JMG));
feenox_check_alloc(neutron_sn.Ai = gsl_matrix_calloc(JMG, JMG));
feenox_check_alloc(neutron_sn.Xi = gsl_matrix_calloc(JMG, JMG));
if (neutron_sn.has_fission) {
feenox_check_alloc(neutron_sn.Fi = gsl_matrix_calloc(JMG, JMG));
}
if (neutron_sn.has_sources) {
feenox_check_alloc(neutron_sn.Si = gsl_vector_calloc(JMG));
}

feenox_check_alloc(neutron_sn.P = gsl_matrix_calloc(MG, JMG));
feenox_check_alloc(neutron_sn.OMEGAB = gsl_matrix_calloc(MG, JMG));
feenox_check_alloc(neutron_sn.DB = gsl_matrix_calloc(MG, JMG));
feenox_check_alloc(neutron_sn.AH = gsl_matrix_calloc(MG, JMG));
if (neutron_sn.has_fission) {
feenox_check_alloc(neutron_sn.XH = gsl_matrix_calloc(MG, JMG));
Expand All @@ -76,12 +79,12 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
material_t *material = feenox_mesh_get_material(e);

// initialize to zero
gsl_matrix_set_zero(neutron_sn.removal);
gsl_matrix_set_zero(neutron_sn.R);
if (neutron_sn.has_fission) {
gsl_matrix_set_zero(neutron_sn.nufission);
gsl_matrix_set_zero(neutron_sn.X);
}
if (neutron_sn.has_sources) {
gsl_vector_set_zero(neutron_sn.src);
gsl_vector_set_zero(neutron_sn.s);
}

// TODO: if XS are uniform then compute only once these cached properties
Expand Down Expand Up @@ -119,28 +122,28 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
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);
gsl_vector_set(neutron_sn.s, 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++) {
for (unsigned int m_prime = 0; m_prime < neutron_sn.directions; m_prime++) {
// scattering
double xi = -neutron_sn.w[n_prime] * neutron_sn.Sigma_s0[g_prime][g].value;
double s = -neutron_sn.w[m_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]);
s -= neutron_sn.w[m_prime] * neutron_sn.Sigma_s1[g_prime][g].value * 3.0 * feenox_mesh_dot(neutron_sn.Omega[m], neutron_sn.Omega[m_prime]);
}
gsl_matrix_set(neutron_sn.removal, diag, dof_index(n_prime,g_prime), xi);
gsl_matrix_set(neutron_sn.R, diag, dof_index(m_prime,g_prime), s);

// 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);
gsl_matrix_set(neutron_sn.X, diag, dof_index(m_prime,g_prime), +neutron_sn.w[m_prime] * neutron_sn.chi[g] * neutron_sn.nu_Sigma_f[g_prime].value);
}
}
}

// absorption
double xi = gsl_matrix_get(neutron_sn.removal, diag, diag);
double xi = gsl_matrix_get(neutron_sn.R, diag, diag);
if (neutron_sn.Sigma_t[g].defined) {
xi += neutron_sn.Sigma_t[g].value;
} else {
Expand All @@ -149,7 +152,7 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
xi += neutron_sn.Sigma_s0[g][g_prime].value;
}
}
gsl_matrix_set(neutron_sn.removal, diag, diag, xi);
gsl_matrix_set(neutron_sn.R, diag, diag, xi);
}
}

Expand All @@ -160,10 +163,10 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
}

// initialize Ki Ai Xi Si <- 0
gsl_matrix_set_zero(neutron_sn.Ki);
gsl_matrix_set_zero(neutron_sn.Li);
gsl_matrix_set_zero(neutron_sn.Ai);
if (neutron_sn.has_fission) {
gsl_matrix_set_zero(neutron_sn.Xi);
gsl_matrix_set_zero(neutron_sn.Fi);
}
if (neutron_sn.has_sources) {
gsl_vector_set_zero(neutron_sn.Si);
Expand All @@ -189,43 +192,43 @@ int feenox_problem_build_volumetric_gauss_point_neutron_sn(element_t *e, unsigne
}

// 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));
feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_sn.D, e->B_G[q], 0.0, neutron_sn.DB));
feenox_call(gsl_blas_dgemm(CblasTrans, CblasNoTrans, e->w[q], neutron_sn.P, neutron_sn.DB, 1.0, neutron_sn.Li));

// petrov-stabilized removal 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(CblasNoTrans, CblasNoTrans, 1.0, neutron_sn.R, 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));
feenox_call(gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, neutron_sn.X, 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.Fi));
}
// 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));
feenox_call(gsl_blas_dgemv(CblasTrans, e->w[q], neutron_sn.P, neutron_sn.s, 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);
gsl_matrix_add(neutron_sn.Li, 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);
gsl_matrix_scale(neutron_sn.Fi, -1.0);
gsl_matrix_add(neutron_sn.Li, neutron_sn.Fi);
} else {
gsl_matrix_add(feenox.pde.Mi, neutron_sn.Xi);
gsl_matrix_add(feenox.pde.Mi, neutron_sn.Fi);
}
}

if (neutron_sn.has_sources) {
gsl_vector_add(feenox.pde.bi, neutron_sn.Si);
}

gsl_matrix_add(feenox.pde.Ki, neutron_sn.Ki);
gsl_matrix_add(feenox.pde.Ki, neutron_sn.Li);
#endif

return FEENOX_OK;
Expand Down
20 changes: 10 additions & 10 deletions src/pdes/neutron_sn/init.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -481,21 +481,21 @@ int feenox_problem_init_runtime_neutron_sn(void) {
}

// 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)));
unsigned int M = neutron_sn.directions;
unsigned int MG = M*G;
feenox_check_alloc(neutron_sn.R = gsl_matrix_calloc(MG, MG));
feenox_check_alloc(neutron_sn.X = gsl_matrix_calloc(MG, MG));
feenox_check_alloc(neutron_sn.s = gsl_vector_calloc(MG));
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++) {
feenox_check_alloc(neutron_sn.D = gsl_matrix_calloc(MG, MG * feenox.pde.dim));
for (unsigned int m = 0; m < M; m++) {
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]);
for (unsigned int d = 0; d < feenox.pde.dim; d++) {
gsl_matrix_set(neutron_sn.D, dof_index(m,g), d*MG + dof_index(m,g), neutron_sn.Omega[m][d]);
}
}
}
Expand Down
23 changes: 15 additions & 8 deletions src/pdes/neutron_sn/neutron_sn.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,13 +65,13 @@ struct neutron_sn_t {


// removal XSs: groups x groups
gsl_matrix *removal;
gsl_matrix *R;

// fission XSs: groups x groups
gsl_matrix *nufission;
gsl_matrix *X;

// independent sources: groups
gsl_vector *src;
gsl_vector *s;

// fission spectrum: groups
double *chi;
Expand All @@ -81,15 +81,22 @@ struct neutron_sn_t {
// elemental matrices
unsigned int n_nodes;

// Petrov-stabilized H_cG
gsl_matrix *P;
gsl_matrix *OMEGA;
gsl_matrix *OMEGAB;

// Direction matrix
gsl_matrix *D;

// intermediate
gsl_matrix *DB;
gsl_matrix *AH;
gsl_matrix *XH;

gsl_matrix *Ki;
gsl_matrix *Ai;
gsl_matrix *Xi;
// elemental matrices (size J*M*G x J*M*G)

gsl_matrix *Li; // leakage
gsl_matrix *Ai; // absorption
gsl_matrix *Fi; // fission
gsl_vector *Si;


Expand Down

0 comments on commit 275f76c

Please sign in to comment.