Initial commit from Polytechnique Montreal

1 files changed, 158 insertions, 0 deletions

diff --git a/Trivac/src/NSSANM1.f90 b/Trivac/src/NSSANM1.f90
new file mode 100755
index 0000000..f16c91a
--- /dev/null
+++ b/Trivac/src/NSSANM1.f90
@@ -0,0 +1,158 @@
+subroutine NSSANM1(nel,ng,nmix,iqfr,qfr,mat,xxx,keff,diff,sigr,chi,sigf,scat,fd,savg)
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Compute the ANM volume fluxes and boundary fluxes and currents using
+! a solution of one- and two-node relations in Cartesian 1D geometry.
+!
+!Copyright:
+! Copyright (C) 2022 Ecole Polytechnique de Montreal
+! This library is free software; you can redistribute it and/or
+! modify it under the terms of the GNU Lesser General Public
+! License as published by the Free Software Foundation; either
+! version 2.1 of the License, or (at your option) any later version
+!
+!Author(s): A. Hebert
+!
+!Parameters: input
+! nel     number of nodes in the nodal calculation.
+! ng      number of energy groups.
+! nmix    number of material mixtures in the nodal calculation.
+! iqfr    node-ordered physical albedo indices.
+! qfr     albedo function information.
+! mat     material mixture index in eacn node.
+! xxx     Cartesian coordinates along the X axis.
+! keff    effective multiplication facctor.
+! diff    diffusion coefficients
+! sigr    removal cross sections.
+! chi     fission spectra.
+! sigf    nu times fission cross section.
+! scat    scattering cross section.
+! fd      discontinuity factors
+! savg    nodal fluxes.
+!
+!Parameters: output
+! savg    boundary fluxes and currents.
+!
+!-----------------------------------------------------------------------
+  !
+  !----
+  !  subroutine arguments
+  !----
+  integer,intent(in) :: nel,ng,nmix,iqfr(6,nel),mat(nel)
+  real,intent(in) :: qfr(6,nel,ng),xxx(nel+1),keff,diff(nmix,ng),sigr(nmix,ng), &
+  & chi(nmix,ng),sigf(nmix,ng),scat(nmix,ng,ng),fd(nmix,2,ng,ng)
+  real, dimension(4*nel+1,ng),intent(inout) :: savg
+  !----
+  ! allocatable arrays
+  !----
+  real, allocatable, dimension(:) :: work1,work2,work4,work5
+  real, allocatable, dimension(:,:) :: A,B,Lambda,work3
+  real(kind=8), allocatable, dimension(:,:,:) :: Lx,Rx
+  !----
+  ! scratch storage allocation
+  !----
+  allocate(A(ng,ng+1),B(ng,ng),Lambda(ng,ng))
+  allocate(work1(ng),work2(ng),work3(ng,ng),work4(ng),work5(ng))
+  allocate(Lx(ng,2*ng,nel),Rx(ng,2*ng,nel))
+  !
+  ! compute nodal coefficients
+  do iel=1,nel
+    ibm=mat(iel)
+    if(ibm == 0) cycle
+    work1(:ng)=diff(ibm,:ng)
+    work2(:ng)=sigr(ibm,:ng)
+    work3(:ng,:ng)=scat(ibm,:ng,:ng)
+    work4(:ng)=chi(ibm,:ng)
+    work5(:ng)=sigf(ibm,:ng)
+    delx=xxx(iel+1)-xxx(iel)
+    call NSSLR1(keff,ng,delx,work1,work2,work3,work4,work5, &
+    & Lx(1,1,iel),Rx(1,1,iel))
+  enddo
+  !----
+  ! compute boundary currents
+  ! left one-node relation
+  !----
+  A(:ng,:ng+1)=0.0
+  if((iqfr(1,1) > 0).or.(iqfr(1,1) == -1)) then
+    ! physical albedo
+    Lambda(:ng,:ng)=0.0
+    do ig=1,ng
+      Lambda(ig,ig)=qfr(1,1,ig)
+    enddo
+    A(:ng,:ng)=real(matmul(Lambda(:ng,:ng),Lx(:ng,ng+1:2*ng,1)),4)
+    B(:ng,:ng)=real(matmul(Lambda(:ng,:ng),Lx(:ng,:ng,1)),4)
+    do ig=1,ng
+      A(ig,ig)=1.0+A(ig,ig)
+    enddo
+    A(:ng,ng+1)=-matmul(B(:ng,:ng),savg(1,:ng))
+  else if(iqfr(1,1) == -2) then
+    ! zero net current
+    do ig=1,ng
+      A(ig,ig)=1.0
+    enddo
+  else if(iqfr(1,1) == -3) then
+    ! zero flux
+    A(:ng,:ng)=real(Lx(:ng,ng+1:2*ng,1),4)
+    A(:ng,ng+1)=real(-matmul(Lx(:ng,:ng,1),savg(1,:ng)),4)
+  else
+    call XABORT('NSSANM1: illegal left boundary condition.')
+  endif
+  call ALSB(ng,1,A,ier,ng)
+  if(ier /= 0) call XABORT('NSSANM1: singular matrix.(1)')
+  savg(3*nel+1,:ng)=A(:ng,ng+1)
+  ! two-node relations
+  do i=2,nel
+    A(:ng,:ng)=real(matmul(fd(mat(i-1),2,:ng,:ng),Rx(:ng,ng+1:2*ng,i-1))- &
+                  & matmul(fd(mat(i),1,:ng,:ng),Lx(:ng,ng+1:2*ng,i)),4)
+    A(:ng,ng+1)=-real(matmul(matmul(fd(mat(i-1),2,:ng,:ng),Rx(:ng,:ng,i-1)),savg(i-1,:ng))- &
+                  & matmul(matmul(fd(mat(i),1,:ng,:ng),Lx(:ng,:ng,i)),savg(i,:ng)),4)
+    call ALSB(ng,1,A,ier,ng)
+    if(ier /= 0) call XABORT('NSSANM1: singular matrix.(2)')
+    savg(3*nel+i,:ng)=A(:ng,ng+1)
+  enddo
+  ! right one-node relation
+  if((iqfr(2,nel) > 0).or.(iqfr(2,nel) == -1)) then
+    ! physical albedo
+    Lambda(:ng,:ng)=0.0
+    do ig=1,ng
+      Lambda(ig,ig)=qfr(2,nel,ig)
+    enddo
+    A(:ng,:ng)=real(matmul(Lambda(:ng,:ng),Rx(:ng,ng+1:2*ng,nel)),4)
+    B(:ng,:ng)=real(matmul(Lambda(:ng,:ng),Rx(:ng,:ng,nel)),4)
+    do ig=1,ng
+      A(ig,ig)=-1.0+A(ig,ig)
+    enddo
+    A(:ng,ng+1)=-matmul(B(:ng,:ng),savg(nel,:ng))
+  else if(iqfr(2,nel) == -2) then
+    ! zero net current
+    do ig=1,ng
+      A(2*nel*ng+ig,2*nel*ng+ig)=1.0
+    enddo
+  else if(iqfr(2,nel) == -3) then
+    ! zero flux
+    A(:ng,:ng)=real(Rx(:ng,ng+1:2*ng,nel),4)
+    A(:ng,ng+1)=real(-matmul(Rx(:ng,:ng,nel),savg(nel,:ng)),4)
+  else
+    call XABORT('NSSANM1: illegal right boundary condition.')
+  endif
+  call ALSB(ng,1,A,ier,ng)
+  if(ier /= 0) call XABORT('NSSANM1: singular matrix.(3)')
+  savg(4*nel+1,:ng)=A(:ng,ng+1)
+  !----
+  ! compute boundary fluxes
+  !----
+  do i=1,nel
+    savg(nel+i,:ng)=real(matmul(Lx(:ng,:ng,i),savg(i,:ng))+ &
+    & matmul(Lx(:ng,ng+1:2*ng,i),savg(3*nel+i,:ng)),4)
+    savg(2*nel+i,:ng)=real(matmul(Rx(:ng,:ng,i),savg(i,:ng))+ &
+    & matmul(Rx(:ng,ng+1:2*ng,i),savg(3*nel+i+1,:ng)),4)
+  enddo
+  !----
+  ! scratch storage deallocation
+  !----
+  deallocate(Rx,Lx)
+  deallocate(work5,work4,work3,work2,work1)
+  deallocate(Lambda,B,A)
+end subroutine NSSANM1