From 5b50993c2ce17aa61a3c0e473418cf4affbf65e8 Mon Sep 17 00:00:00 2001
From: Wiarda <wiardada@ornl.gov>
Date: Thu, 16 Sep 2021 10:37:16 -0400
Subject: [PATCH] Make the calcCross function a procedure on
CrossSectionCalculator. This way we can call it as a normal virtual function.
Let CrossSectionCalcDriver only have one CrossSectionCalculator
Unfortunately that needs to be a pointer to the
R-Matrix formalism that is actually used.
---
sammy/src/clq/ArtificalCross_M.f90 | 9 +-
sammy/src/mlb/mmlb2.f90 | 46 ++--
sammy/src/mlb/mmlb3.f90 | 6 +-
sammy/src/mlb/mmlb4.f90 | 4 +-
sammy/src/the/CrossSectionCalcDriver_M.f90 | 249 +++++++--------------
sammy/src/the/CrossSectionCalculator_M.f90 | 20 ++
6 files changed, 145 insertions(+), 189 deletions(-)
diff --git a/sammy/src/clq/ArtificalCross_M.f90 b/sammy/src/clq/ArtificalCross_M.f90
index c9ed83730..4de7e5290 100644
--- a/sammy/src/clq/ArtificalCross_M.f90
+++ b/sammy/src/clq/ArtificalCross_M.f90
@@ -22,7 +22,7 @@ module ArtificalCross_M
procedure, pass(this) :: initialize => ArtificalCross_initialize
procedure, pass(this) :: destroy => ArtificalCross_destroy
procedure, pass(this) :: setRange => ArtificalCross_setRange
- procedure, pass(this) :: calculate => ArtificalCross_calculate
+ procedure, pass(this) :: calcCross => ArtificalCross_calcCross
procedure, pass(this) :: print => ArtificalCross_print
end type ArtificalCross
@@ -156,9 +156,10 @@ subroutine ArtificalCross_print(this)
END IF
end subroutine
-subroutine ArtificalCross_calculate(this, ener)
- class(ArtificalCross) :: this
- real(kind=8)::ener
+subroutine ArtificalCross_calcCross(this, ener, Ipoten)
+ class(ArtificalCross) :: this
+ real(kind=8)::ener
+ integer::ipoten
real(kind=8)::Su, Squ, Sigxxx, X
diff --git a/sammy/src/mlb/mmlb2.f90 b/sammy/src/mlb/mmlb2.f90
index de615e707..44243fded 100755
--- a/sammy/src/mlb/mmlb2.f90
+++ b/sammy/src/mlb/mmlb2.f90
@@ -1,39 +1,51 @@
-module mmlb2_m
+module MlbwCrossCalcImpl_m
use MlbwCrossCalc_M
+ use CrossSectionCalculator_M
use mmlb4_m
implicit none
+
+ type, extends(MlbwCrossCalc) :: MlbwCrossCalcImpl
+ contains
+ procedure, pass(this) :: calcCross => MlbwCrossCalcImpl_calcCross
+ end type
+
contains
- subroutine MlbwCrossCalc_calculateCross(calculator, ener, Ipoten)
- class(MlbwCrossCalc) :: calculator
+ subroutine MlbwCrossCalcImpl_calcCross(this, ener, Ipoten)
+ class(MlbwCrossCalcImpl) ::this
real(kind=8)::ener
integer::ipoten
real(kind=8)::val
integer::i, itot
+ logical(C_BOOL)::accu
- itot = calculator%covariance%getNumTotalParam()
- call calculator%crossData%reserveColumns(calculator%row, itot + 1)
+ call CrossSectionCalculator_calcCross(this, ener, Ipoten)
if (ener.eq.0.0d0) return
- calculator%ener = ener
- IF (calculator%ener.LT.0.0d0) calculator%ener = -calculator%ener
- calculator%enerSq = dSQRT(calculator%ener)
+ IF (this%ener.LT.0.0d0) this%ener = -this%ener
- CALL Abpart_Mlb (calculator)
- CALL Parsh_Mlb ( calculator, Ipoten)
+ CALL Abpart_Mlb (this)
+ CALL Parsh_Mlb (this, Ipoten)
- if (ener.lt.0.0d0) then
+ if (ener.lt.0.0d0) then
+ itot = this%covariance%getNumTotalParam() + 1
+ if ( .not.this%wantDerivs) itot = 1
+
+ accu = .false.
+ call this%crossData%setAccumulate(accu)
do i = 0, itot
- val = calculator%crossData%getData(calculator%row, i, 1)
+ val = this%crossData%getData(this%row, i, 1)
if (val.ne.0.0d0) then
- ! addData accumulates. So to change the sign,
- ! subtract twice the original value
- val = -2.0d0*val
- call calculator%crossData%addData(calculator%row, i, 1, val)
+ val = -val
+ call this%crossData%addData(this%row, i, 1, val)
end if
end do
+ accu = .true.
+ call this%crossData%setAccumulate(accu)
end if
+
+ this%ener = ener
end subroutine
-end module mmlb2_m
+end module MlbwCrossCalcImpl_m
diff --git a/sammy/src/mlb/mmlb3.f90 b/sammy/src/mlb/mmlb3.f90
index 876c996c3..5f5e2c743 100644
--- a/sammy/src/mlb/mmlb3.f90
+++ b/sammy/src/mlb/mmlb3.f90
@@ -15,7 +15,7 @@ contains
!
use constn_common_m, only : Pi100
!
- type(MlbwCrossCalc)::calculator
+ class(MlbwCrossCalc)::calculator
real(kind=8),intent(in)::Goj
integer, intent(in)::Nent, Next, If_Zke, Iso, igr, ipar
@@ -127,7 +127,7 @@ contains
!
real(kind=8),intent(in)::Goj
integer,intent(in)::Nent, Next, If_Zke, igr, ipar
- type(MlbwCrossCalc)::calculator
+ class(MlbwCrossCalc)::calculator
real(kind=8)::val, Zero, Sum, Sum1, Sum2
real(kind=8)::A,B,C
@@ -196,7 +196,7 @@ contains
use constn_common_m, only : Pi100
use SammySpinGroupInfo_M
!
- type(MlbwCrossCalc)::calculator
+ class(MlbwCrossCalc)::calculator
type(SammySpinGroupInfo)::spinInfo
integer::imin, ipar
diff --git a/sammy/src/mlb/mmlb4.f90 b/sammy/src/mlb/mmlb4.f90
index 9606498f5..2285c3fbe 100644
--- a/sammy/src/mlb/mmlb4.f90
+++ b/sammy/src/mlb/mmlb4.f90
@@ -21,7 +21,7 @@ module mmlb4_m
use RMatResonanceParam_M
use xxx3
!
- type(MlbwCrossCalc)::calculator
+ class(MlbwCrossCalc)::calculator
type(SammyResonanceInfo)::resInfo
type(SammySpinGroupInfo)::spinInfo
type(RMatResonance)::resonance
@@ -233,7 +233,7 @@ module mmlb4_m
use xxx6
!
!
- type(MlbwCrossCalc)::calculator
+ class(MlbwCrossCalc)::calculator
integer::Ipoten
type(SammySpinGroupInfo)::spinInfo
diff --git a/sammy/src/the/CrossSectionCalcDriver_M.f90 b/sammy/src/the/CrossSectionCalcDriver_M.f90
index 73c55515e..2b08c3700 100644
--- a/sammy/src/the/CrossSectionCalcDriver_M.f90
+++ b/sammy/src/the/CrossSectionCalcDriver_M.f90
@@ -1,6 +1,6 @@
module CrossSectionCalcDriver_M
use CrossSectionCalculator_M
- use MlbwCrossCalc_M
+ use MlbwCrossCalcImpl_m
use CroCrossCalc_M
use XctCrossCalc_M
use ArtificalCross_M
@@ -12,23 +12,18 @@ module CrossSectionCalcDriver_M
type CrossSectionCalcDriver
- ! what kind of cross section do we want
- type(MlbwCrossCalc)::mlbw ! multi-level (or single level) Breit-Wigner
- logical::mlbwInit=.false.
-
- type(ArtificalCross)::artifical ! artificial cross section data for testing
- logical::artificalInit=.false.
-
- type(CroCrossCalc)::cro ! old style Reich-Moore (no multi-isotope no angular)
- logical::croInit=.false.
-
- type(XctCrossCalc)::xct ! new style Reich-Moore
- logical::xctInit=.false.
-
! arrays to take the cross section and derivative data
type(DerivativeHandler)::calcData
type(DerivativeHandler)::calcDataSelf ! used if self-indication experiments are calculated
+ class(CrossSectionCalculator),pointer::calculator
+ logical::calculatorInit
+
+ type(ArtificalCross),pointer::artificalCalc => null()
+ type(MlbwCrossCalcImpl),pointer,private::mlbwCalc => null()
+ type(CroCrossCalc),pointer,private::croCalc => null()
+ type(XctCrossCalc),pointer,private::xctCalc => null()
+
logical::wantCrossPerSpin = .false. ! do we want cross section per spin group
type(GridData)::crossPerSpin ! if used: energy, group 1, group 1 + group 2, group 1 + group 2 + group 3, ... repeated for each Nnnsig value )
@@ -51,7 +46,7 @@ module CrossSectionCalcDriver_M
procedure, pass(this) :: setWantCrossPerSpin => CrossSectionCalcDriver_setWantCrossPerSpin ! indicate whether we want cross section per spin
end type CrossSectionCalcDriver
- private setup_desiredCross, setUpCommonData
+ private setup_desiredCross
contains
subroutine CrossSectionCalcDriver_initialize(this, calcData, calcDataSelf, kwcoul, Kssmsc, Kcros)
@@ -87,61 +82,66 @@ module CrossSectionCalcDriver_M
IF (this%Kcros.EQ.11) needAngular = .true.
first = .false.
- if (Kkkclq.NE.0) THEN
- if (.not.this%artificalInit) then
- call this%artifical%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
- this%artifical%Kkclqx = Kkclqx
- this%artifical%Kkkclq = Kkkclq
-
- this%artificalInit = .true.
- first = .true.
- end if
- call setUpCommonData(this, this%artifical, first, this%Kcros, covData, resParData, radFitFlags, wantDeriv)
- ELSE IF (Krmatx.EQ.1 .OR. Krmatx.EQ.-1) THEN
- if (.not.this%mlbwInit) then
- call this%mlbw%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
- this%mlbw%wantMlbw = .false.
- if (Krmatx.eq.1) then
- this%mlbw%wantMlbw = .true.
- else if (Krmatx.eq.-1) then
- this%mlbw%wantMlbw = .false.
- else
- STOP '[STOP in mlb/mmlb0.f expect MLBW or SLBW]'
- end if
-
- this%mlbwInit = .true.
- first = .true.
- end if
- call setUpCommonData(this, this%mlbw, first, this%Kcros, covData, resParData, radFitFlags, wantDeriv)
- ELSE IF (Krmatx.EQ.0) THEN
- if (.not.this%croInit) then
- call this%cro%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
-
- this%croInit = .true.
- first = .true.
- end if
- call setUpCommonData(this, this%cro, first, this%Kcros, covData, resParData, radFitFlags, wantDeriv)
- ELSE IF (Krmatx.EQ.2) THEN
- if (.not.this%xctInit) then
- nisoOur = niso
- if (needAngular) then
- call info%initialize(resParData, covData, radFitFlags)
- iq = info%getNumUniqEchanIncludedChannels()
- if( iq.gt.nisoOur) nisoOur = iq
- call info%destroy()
- end if
-
- call this%xct%initialize(resParData, covData, radFitFlags, nisoOur, needAngular, Itzero, Ilzero)
-
- this%xctInit = .true.
- first = .true.
+ nisoOur = 1
+ if (.not.this%calculatorInit) then
+ if (Kkkclq.NE.0) THEN
+ allocate(this%artificalCalc)
+ call this%artificalCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
+ this%artificalCalc%Kkclqx = Kkclqx
+ this%artificalCalc%Kkkclq = Kkkclq
+ this%calculator => this%artificalCalc
+ ELSE IF (Krmatx.EQ.1 .OR. Krmatx.EQ.-1) THEN
+ allocate(this%mlbwCalc)
+ call this%mlbwCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
+ if (Krmatx.eq.1) then
+ this%mlbwCalc%wantMlbw = .true.
+ else if (Krmatx.eq.-1) then
+ this%mlbwCalc%wantMlbw = .false.
+ else
+ STOP '[STOP in mlb/mmlb0.f expect MLBW or SLBW]'
+ end if
+ this%calculator => this%mlbwCalc
+ ELSE IF (Krmatx.EQ.0) THEN
+ allocate(this%croCalc)
+ call this%croCalc%initialize(resParData, covData, radFitFlags, 1, needAngular, Itzero, Ilzero)
+ this%calculator => this%croCalc
+ ELSE IF (Krmatx.EQ.2) THEN
+ nisoOur = niso
+ if (needAngular) then
+ call info%initialize(resParData, covData, radFitFlags)
+ iq = info%getNumUniqEchanIncludedChannels()
+ if( iq.gt.nisoOur) nisoOur = iq
+ call info%destroy()
+ end if
+ allocate(this%xctCalc)
+ call this%xctCalc%initialize(resParData, covData, radFitFlags, nisoOur, needAngular, Itzero, Ilzero)
+ this%calculator => this%xctCalc
+ ELSE
+ write(0,*)" Invalid selection for cross section calculator"
+ stop
end if
- call setUpCommonData(this, this%xct, first, this%Kcros, covData, resParData, radFitFlags, wantDeriv)
+ this%calculatorInit = .true.
+ first = .true.
end if
+
+ this%calculator%crossData%instance_ptr = this%calcData%instance_ptr
+ if (first) then
+ call this%calculator%crossData%setUpList(resParData, radFitFlags, nisoOur)
+ else
+ ! unfortunately they sometimes get redefined
+ this%calculator%resData%instance_ptr = resParData%instance_ptr
+ this%calculator%covariance%instance_ptr = covData%instance_ptr
+ this%calculator%radiusData%instance_ptr = radFitFlags%instance_ptr
+ end if
+
+ this%calculator%reactType = this%Kcros
+ call setup_desiredCross(this, this%calculator)
+ call this%calculator%setWantDeriv(wantDeriv)
+ call this%calculator%setUpDerivativeList()
+ call this%calculator%setEnergyIndependent(this%Kcros, Twomhb, this%kwcoul, Etac)
end subroutine
subroutine CrossSectionCalcDriver_calcCross(this, row, ener, Ipoten)
- use mmlb2_m
class(CrossSectionCalcDriver)::this
real(kind=8)::ener
integer::ipoten, row
@@ -153,21 +153,9 @@ module CrossSectionCalcDriver_M
call this%calcData%setAccumulate(empty)
- if (this%mlbwInit) then
- this%mlbw%row = row
- ! the specific call is here instead
- ! a call function for the this%mlbw, because
- ! MlbwCrossCalc_calculateCross relies
- ! on data in this%mlbw
- ! this avoids a circular dependence
- call MlbwCrossCalc_calculateCross(this%mlbw, ener, Ipoten)
- else if (this%artificalInit) then
- this%artifical%row = row
- call this%artifical%calculate(ener)
- else if (this%croInit) then
- this%cro%row = row
- else if (this%xctInit) then
- this%xct%row = row
+ if (this%calculatorInit) then
+ this%calculator%row = row
+ call this%calculator%calcCross( ener, Ipoten)
else
write(0,*)" Cross section calculator not initialized"
stop
@@ -179,33 +167,6 @@ module CrossSectionCalcDriver_M
end subroutine
- subroutine setUpCommonData(this, calc, first, Kcros, covData, resParData, radFitFlags, wantDeriv)
- class(CrossSectionCalcDriver)::this
- class(CrossSectionCalculator)::calc
- logical::first, wantDeriv
- integer::Kcros
- type(SammyRMatrixParameters)::resParData
- type(CovarianceData)::covData
- type(AdjustedRadiusData)::radFitFlags
-
- calc%crossData%instance_ptr = this%calcData%instance_ptr
- if (first) then
- call calc%crossData%setUpList(resParData, radFitFlags, calc%numIso)
- else
- ! unfortunately they sometimes get redefined
- calc%resData%instance_ptr = resParData%instance_ptr
- calc%covariance%instance_ptr = covData%instance_ptr
- calc%radiusData%instance_ptr = radFitFlags%instance_ptr
- end if
-
- calc%reactType = Kcros
- call setup_desiredCross(this, calc)
- call calc%setWantDeriv(wantDeriv)
- call calc%setUpDerivativeList()
- call calc%setEnergyIndependent(Kcros, Twomhb, this%kwcoul, Etac)
- end subroutine setUpCommonData
-
-
subroutine setup_desiredCross(this, calc)
use fixedi_m, only : Kaptur
use ifwrit_m, only : Kaverg, Kfinit, Krecon
@@ -292,33 +253,22 @@ module CrossSectionCalcDriver_M
subroutine CrossSectionCalcDriver_destroy(this)
class(CrossSectionCalcDriver)::this
- if (this%mlbwInit) call this%mlbw%destroy()
- this%mlbwInit = .false.
- if (this%artificalInit) call this%artifical%destroy()
- this%artificalInit = .false.
- if (this%croInit) call this%cro%destroy()
- this%croInit = .false.
- if (this%xctInit) call this%xct%destroy()
- this%xctInit = .false.
-
- if (this%wantCrossPerSpin) then
- call this%crossPerSpin%destroy()
- this%wantCrossPerSpin = .false.
+ if (this%calculatorInit) then
+ call this%calculator%destroy()
end if
+ this%calculatorInit = .false.
+ if (associated(this%artificalCalc)) deallocate(this%artificalCalc)
+ if (associated(this%mlbwCalc)) deallocate(this%mlbwCalc)
+ if (associated(this%croCalc)) deallocate(this%croCalc)
+ if (associated(this%xctCalc)) deallocate(this%xctCalc)
end subroutine
subroutine CrossSectionCalcDriver_getAbundanceFromSpinGroup(this, setAbund)
class(CrossSectionCalcDriver)::this
logical::setAbund
- if (this%mlbwInit) then
- this%mlbw%getAbundanceFromSpinGroup = setAbund
- else if (this%artificalInit) then
- this%artifical%getAbundanceFromSpinGroup = setAbund
- else if (this%croInit) then
- this%cro%getAbundanceFromSpinGroup = setAbund
- else if (this%xctInit) then
- this%xct%getAbundanceFromSpinGroup = setAbund
+ if (this%calculatorInit) then
+ this%calculator%getAbundanceFromSpinGroup = setAbund
end if
end subroutine
@@ -327,14 +277,8 @@ module CrossSectionCalcDriver_M
real(kind=8)::emin,emax
integer::num
- if (this%mlbwInit) then
- call this%mlbw%setRange(emin, emax, num)
- else if (this%artificalInit) then
- call this%artifical%setRange(emin, emax, num)
- else if (this%croInit) then
- call this%cro%setRange(emin, emax, num)
- else if (this%xctInit) then
- call this%xct%setRange(emin, emax, num)
+ if (this%calculatorInit) then
+ call this%calculator%setRange(emin, emax, num)
end if
end subroutine
@@ -342,14 +286,8 @@ module CrossSectionCalcDriver_M
class(CrossSectionCalcDriver) :: this
value = .false.
- if (this%mlbwInit) then
- value = this%mlbw%wantDerivs
- else if (this%artificalInit) then
- value = this%artifical%wantDerivs
- else if (this%croInit) then
- value = this%cro%wantDerivs
- else if (this%xctInit) then
- value = this%xct%wantDerivs
+ if (this%calculatorInit) then
+ value = this%calculator%wantDerivs
end if
end function
@@ -358,14 +296,8 @@ module CrossSectionCalcDriver_M
value = 0
- if (this%mlbwInit) then
- value = this%mlbw%covariance%getNumTotalParam()
- else if (this%artificalInit) then
- value = this%artifical%covariance%getNumTotalParam()
- else if (this%croInit) then
- value = this%cro%covariance%getNumTotalParam()
- else if (this%xctInit) then
- value = this%xct%covariance%getNumTotalParam()
+ if (this%calculatorInit) then
+ value = this%calculator%covariance%getNumTotalParam()
end if
end function
@@ -382,18 +314,9 @@ module CrossSectionCalcDriver_M
call this%crossPerSpin%initialize()
instance_ptr =this%crossPerSpin%instance_ptr
end if
- if (this%mlbwInit) then
- this%mlbw%wantCrossPerSpin = want
- this%mlbw%crossPerSpin%instance_ptr = instance_ptr
- else if (this%artificalInit) then
- this%artifical%wantCrossPerSpin = want
- this%artifical%crossPerSpin%instance_ptr = instance_ptr
- else if (this%croInit) then
- this%cro%wantCrossPerSpin = want
- this%cro%crossPerSpin%instance_ptr = instance_ptr
- else if (this%xctInit) then
- this%xct%wantCrossPerSpin = want
- this%xct%crossPerSpin%instance_ptr = instance_ptr
+ if (this%calculatorInit) then
+ this%calculator%wantCrossPerSpin = want
+ this%calculator%crossPerSpin%instance_ptr = instance_ptr
end if
end subroutine
end module CrossSectionCalcDriver_M
diff --git a/sammy/src/the/CrossSectionCalculator_M.f90 b/sammy/src/the/CrossSectionCalculator_M.f90
index 5ccbf0536..fd4aea11a 100644
--- a/sammy/src/the/CrossSectionCalculator_M.f90
+++ b/sammy/src/the/CrossSectionCalculator_M.f90
@@ -69,6 +69,7 @@ module CrossSectionCalculator_M
procedure, pass(this) :: setRange => CrossSectionCalculator_setRange ! set the range and reserve enough space in crossData (the latter for efficiency only)
procedure, pass(this) :: initialize => CrossSectionCalculator_initialize
procedure, pass(this) :: destroy => CrossSectionCalculator_destroy
+ procedure, pass(this) :: calcCross => CrossSectionCalculator_calcCross
end type CrossSectionCalculator
contains
@@ -150,6 +151,25 @@ subroutine CrossSectionCalculator_setWantDeriv(this, wantDerivs)
end if
end subroutine
+subroutine CrossSectionCalculator_calcCross(this, ener, Ipoten)
+ class(CrossSectionCalculator) :: this
+ real(kind=8)::ener
+ integer::ipoten
+
+ integer::itot
+
+
+ itot = this%covariance%getNumTotalParam() + 1
+ if ( .not.this%wantDerivs) itot = 1
+ call this%crossData%reserveColumns(this%row, itot)
+
+ this%ener = ener
+
+ this%enerSq = ener
+ IF (this%ener.LT.0.0d0) this%enerSq = -this%enerSq
+ if (this%enerSq.gt.0.0d0) this%enerSq = dSQRT(this%enerSq)
+end subroutine
+
!
! For a given spin group, determine the number of varied parameters
! (including their resonance and channel.
--
GitLab