Initial commit from Polytechnique Montreal

1 files changed, 219 insertions, 0 deletions

diff --git a/Trivac/src/NSSMXYZ.f90 b/Trivac/src/NSSMXYZ.f90
new file mode 100755
index 0000000..92a9db7
--- /dev/null
+++ b/Trivac/src/NSSMXYZ.f90
@@ -0,0 +1,219 @@
+subroutine NSSMXYZ(ll4f,ndim,nx,ny,nz,nmix,mat,xx,yy,zz,idl,vol,iqfr,qfr, &
+& diff,drift,sigt,mux,muy,muz,imax,imay,imaz,ipy,ipz,a11x,a11y,a11z)
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! Assembly of system matrices for coarse mesh finite differences with
+! nodal correction.
+!
+!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
+! ll4f    total number of averaged flux unknown per energy group.
+! ndim    number of dimensions (1, 2, or 3).
+! nx      number of nodes in the X direction.
+! ny      number of nodes in the Y direction.
+! nz      number of nodes in the Z direction.
+! nmix    number of mixtures.
+! mat     node mixtures.
+! xx      node widths in the X direction.
+! yy      node widths in the Y direction.
+! zz      node widths in the Z direction.
+! idl     position of averaged fluxes in unknown vector.
+! vol     node volumes.
+! iqfr    boundary conditions.
+! qfr     albedo functions.
+! diff    diffusion coefficients.
+! drift   drift coefficients.
+! sigt    removal macroscopic cross section.
+! mux     X-oriented compressed storage mode indices.
+! muy     Y-oriented compressed storage mode indices.
+! muz     Z-oriented compressed storage mode indices.
+! imax    X-oriented position of each first non-zero column element.
+! imay    Y-oriented position of each first non-zero column element.
+! imaz    Z-oriented position of each first non-zero column element.
+! ipy     Y-oriented permutation matrices.
+! ipz     Z-oriented permutation matrices.
+!
+!Parameters: output
+! a11x    X-directed  matrices corresponding to the divergence (i.e.
+!         leakage) and removal terms. Dimensionned to imax(ll4f).
+! a11y    Y-directed  matrices corresponding to the divergence (i.e.
+!         leakage) and removal terms. Dimensionned to imay(ll4f).
+! a11z    Z-directed  matrices corresponding to the divergence (i.e.
+!         leakage) and removal terms. Dimensionned to imaz(ll4f).
+!
+!-----------------------------------------------------------------------
+!
+  !----
+  !  subroutine arguments
+  !----
+  integer,intent(in) :: ll4f,ndim,nx,ny,nz,nmix,mat(nx,ny,nz),idl(nx,ny,nz), &
+  & iqfr(6,nx,ny,nz),mux(ll4f),muy(ll4f),muz(ll4f),imax(ll4f),imay(ll4f),imaz(ll4f), &
+  & ipy(ll4f),ipz(ll4f)
+  real,intent(in) :: xx(nx,ny,nz),yy(nx,ny,nz),zz(nx,ny,nz),vol(nx,ny,nz), &
+  & qfr(6,nx,ny,nz),diff(nmix),drift(6,nx,ny,nz),sigt(nmix)
+  real,intent(out) :: a11x(*),a11y(*),a11z(*)
+  !----
+  !  local variables
+  !----
+  real :: coef(6),codr(6)
+  !
+  a11x(:imax(ll4f))=0.0
+  if(ndim > 1) a11y(:imay(ll4f))=0.0
+  if(ndim == 3) a11z(:imaz(ll4f))=0.0
+  do k=1,nz
+    do j=1,ny
+      do i=1,nx
+        ibm=mat(i,j,k)
+        if(ibm <= 0) cycle
+        kel=idl(i,j,k)
+        if(kel == 0) cycle
+        vol0=vol(i,j,k)
+        call NSSCO(nx,ny,nz,nmix,i,j,k,mat,xx,yy,zz,diff,iqfr(1,i,j,k),qfr(1,i,j,k),coef)
+        coef(1:2)=coef(1:2)*vol0/xx(i,j,k)
+        coef(3:4)=coef(3:4)*vol0/yy(i,j,k)
+        coef(5:6)=coef(5:6)*vol0/zz(i,j,k)
+        codr(1:2)=drift(1:2,i,j,k)*vol0/xx(i,j,k)
+        codr(3:4)=drift(3:4,i,j,k)*vol0/yy(i,j,k)
+        codr(5:6)=drift(5:6,i,j,k)*vol0/zz(i,j,k)
+        !
+        ! x-directed couplings
+        kel2=0
+        kk1=iqfr(1,i,j,k)
+        if(kk1 == -4) then
+          kel2=idl(nx,j,k)
+        else if(kk1 == 0) then
+          kel2=idl(i-1,j,k)
+        endif
+        if(kel2 /= 0) then
+          if(kel2 <= kel) then
+            key=mux(kel)-kel+kel2
+            a11x(key)=a11x(key)-coef(1)+codr(1)
+          else
+            key=mux(kel2)+kel2-kel
+            a11x(key)=a11x(key)-coef(1)+codr(1)
+          endif
+        endif
+        kel2=0
+        kk2=iqfr(2,i,j,k)
+        if(kk2 == -4) then
+          kel2=idl(1,j,k)
+        else if(kk2 == 0) then
+          kel2=idl(i+1,j,k)
+        endif
+        if(kel2 /= 0) then
+          if(kel2 <= kel) then
+            key=mux(kel)-kel+kel2
+            a11x(key)=a11x(key)-coef(2)-codr(2)
+          else
+            key=mux(kel2)+kel2-kel
+            a11x(key)=a11x(key)-coef(2)-codr(2)
+          endif
+        endif
+        key0=mux(kel)
+        a11x(key0)=a11x(key0)+coef(1)+codr(1)+coef(2)-codr(2)
+        a11x(key0)=a11x(key0)+coef(3)+codr(3)+coef(4)-codr(4)
+        a11x(key0)=a11x(key0)+coef(5)+codr(5)+coef(6)-codr(6)
+        a11x(key0)=a11x(key0)+sigt(ibm)*vol0
+        !
+        if(ndim > 1) then
+          ! y-directed couplings
+          kel2=0
+          kk3=iqfr(3,i,j,k)
+          if(kk3 == -4) then
+            kel2=idl(i,ny,k)
+          else if(kk3 == 0) then
+            kel2=idl(i,j-1,k)
+          endif
+          ind1=ipy(kel)
+          if(kel2 /= 0) then
+            ind2=ipy(kel2)
+            if(kel2 <= kel) then
+              key=muy(ind1)-ind1+ind2
+              a11y(key)=a11y(key)-coef(3)+codr(3)
+            else
+              key=muy(ind2)+ind2-ind1
+              a11y(key)=a11y(key)-coef(3)+codr(3)
+            endif
+          endif
+          kel2=0
+          kk4=iqfr(4,i,j,k)
+          if(kk4 == -4) then
+            kel2=idl(i,1,k)
+          else if(kk4 == 0) then
+            kel2=idl(i,j+1,k)
+          endif
+          if(kel2 /= 0) then
+            ind2=ipy(kel2)
+            if(kel2 <= kel) then
+              key=muy(ind1)-ind1+ind2
+              a11y(key)=a11y(key)-coef(4)-codr(4)
+            else
+              key=muy(ind2)+ind2-ind1
+              a11y(key)=a11y(key)-coef(4)-codr(4)
+            endif
+          endif
+          key0=muy(ind1)
+          a11y(key0)=a11y(key0)+coef(1)+codr(1)+coef(2)-codr(2)
+          a11y(key0)=a11y(key0)+coef(3)+codr(3)+coef(4)-codr(4)
+          a11y(key0)=a11y(key0)+coef(5)+codr(5)+coef(6)-codr(6)
+          a11y(key0)=a11y(key0)+sigt(ibm)*vol0
+        endif
+        !
+        if(ndim > 2) then
+          ! z-directed couplings
+          kel2=0
+          kk5=iqfr(5,i,j,k)
+          if(kk5 == -4) then
+            kel2=idl(i,j,nz)
+          else if(kk5 == 0) then
+            kel2=idl(i,j,k-1)
+          endif
+          ind1=ipz(kel)
+          if(kel2 /= 0) then
+            ind2=ipz(kel2)
+            if(kel2 <= kel) then
+              key=muz(ind1)-ind1+ind2
+              a11z(key)=a11z(key)-coef(5)+codr(5)
+            else
+              key=muz(ind2)+ind2-ind1
+              a11z(key)=a11z(key)-coef(5)+codr(5)
+            endif
+          endif
+          kel2=0
+          kk6=iqfr(6,i,j,k)
+          if(kk6 == -4) then
+            kel2=idl(i,j,1)
+          else if(kk6 == 0) then
+            kel2=idl(i,j,k+1)
+          endif
+          if(kel2 /= 0) then
+            ind2=ipz(kel2)
+            if(kel2 <= kel) then
+              key=muz(ind1)-ind1+ind2
+              a11z(key)=a11z(key)-coef(6)-codr(6)
+            else
+              key=muz(ind2)+ind2-ind1
+              a11z(key)=a11z(key)-coef(6)-codr(6)
+            endif
+          endif
+          key0=muz(ind1)
+          a11z(key0)=a11z(key0)+coef(1)+codr(1)+coef(2)-codr(2)
+          a11z(key0)=a11z(key0)+coef(3)+codr(3)+coef(4)-codr(4)
+          a11z(key0)=a11z(key0)+coef(5)+codr(5)+coef(6)-codr(6)
+          a11z(key0)=a11z(key0)+sigt(ibm)*vol0
+        endif
+      enddo
+    enddo
+  enddo
+  return
+end subroutine NSSMXYZ