Skip to content

Commit

Permalink
Updates in SDAlgorithms
Browse files Browse the repository at this point in the history 
- adding HHT alpha and Collocation methods
- working on the issue #442
  • Loading branch information
@shishiousan
shishiousan committed Feb 3, 2024
1 parent 4d16ba3 commit 50e0067
Showing 1 changed file with 171 additions and 1 deletion.
172 changes: 171 additions & 1 deletion src/modules/AbstractKernel/src/SDAlgorithms.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,13 +104,19 @@ MODULE SDAlgorithms
!! acc(4) coefficient of solution
LOGICAL(LGT) :: acc_zero(4) = .TRUE.

REAL(DFP) :: alpha = 0.0_DFP
!! only in case of HHT-alpha method it can be less than 0.0

CONTAINS
PRIVATE

PROCEDURE, PUBLIC, PASS(obj) :: DEALLOCATE => obj_Deallocate
PROCEDURE, PUBLIC, PASS(obj) :: Display => obj_Display
PROCEDURE, PUBLIC, PASS(obj) :: MakeZeros => obj_MakeZeros
PROCEDURE, PUBLIC, PASS(obj) :: NewmarkBeta => obj_NewmarkBeta
PROCEDURE, PUBLIC, PASS(obj) :: HHTAlpha => obj_HHTAlpha
PROCEDURE, PUBLIC, PASS(obj) :: Collocation => obj_Collocation
! PROCEDURE, PUBLIC, PASS(obj) :: Houbolt => obj_Houbolt
PROCEDURE, PASS(obj) :: ImportFromToml1 => obj_ImportFromToml1
PROCEDURE, PASS(obj) :: ImportFromToml2 => obj_ImportFromToml2
GENERIC, PUBLIC :: ImportFromToml => ImportFromToml1, ImportFromToml2
Expand Down Expand Up @@ -173,6 +179,116 @@ SUBROUTINE obj_NewmarkBeta(obj, beta, gamma)

END SUBROUTINE obj_NewmarkBeta

!----------------------------------------------------------------------------
! HHTAlpha
!----------------------------------------------------------------------------

SUBROUTINE obj_HHTAlpha(obj, alpha, beta, gamma)
CLASS(SDAlgoParam_), INTENT(INOUT) :: obj
REAL(DFP), OPTIONAL, INTENT(IN) :: alpha
!! Default -0.30
REAL(DFP), OPTIONAL, INTENT(IN) :: beta
!! Default is (1-alpha)**2/4
REAL(DFP), OPTIONAL, INTENT(IN) :: gamma
!! Default is (1-2alpha)/2

! internal varibales
REAL(DFP) :: alpha0, beta0, gamma0, areal

CALL obj%DEALLOCATE()

alpha0 = Input(default=-0.30_DFP, option=alpha)
areal = (1.0_DFP - alpha0)**2 * 0.25_DFP
beta0 = Input(default=areal, option=beta)
areal = (1.0_DFP - 2.0_DFP * alpha0) * 0.50_DFP
gamma0 = Input(default=areal, option=gamma)

obj%alpha = alpha0

CALL obj%NewmarkBeta(beta=beta0, gamma=gamma0)

obj%tanmat(2) = obj%tanmat(2) * (1.0_DFP + alpha0)
obj%tanmat(3) = obj%tanmat(3) * (1.0_DFP + alpha0)

obj%rhs_u1(2) = obj%rhs_u1(2) * (1.0_DFP + alpha0)
obj%rhs_u1(3) = alpha0 * beta0

obj%rhs_v1(2) = obj%rhs_v1(2) * (1.0_DFP + alpha0) &
& + alpha0 * beta0

obj%rhs_a1(2) = obj%rhs_a1(2) * (1.0_DFP + alpha0)

CALL obj%MakeZeros()

END SUBROUTINE obj_HHTAlpha

!----------------------------------------------------------------------------
! Collocation
!----------------------------------------------------------------------------

SUBROUTINE obj_Collocation(obj, beta, gamma, theta)
CLASS(SDAlgoParam_), INTENT(INOUT) :: obj
REAL(DFP), OPTIONAL, INTENT(IN) :: beta
!! Default is 1/6
REAL(DFP), OPTIONAL, INTENT(IN) :: gamma
!! Default is 0.50
REAL(DFP), OPTIONAL, INTENT(IN) :: theta
!! Default is 1.4 (Wilson theta)

! internal varibales
REAL(DFP) :: theta0, beta0, gamma0, areal

CALL obj%DEALLOCATE()

beta0 = Input(default=1.0_DFP / 6.0_DFP, option=beta)
gamma0 = Input(default=0.50_DFP, option=gamma)
theta0 = Input(default=1.4_DFP, option=theta)

CALL obj%NewmarkBeta(beta=beta0, gamma=gamma0)

obj%tanmat(2) = obj%tanmat(2) * theta0
obj%tanmat(3) = obj%tanmat(3) * theta0**2

obj%rhs_f2 = obj%rhs_f2 * theta0

obj%rhs_u1(2) = obj%rhs_u1(2) * theta0

obj%rhs_v1(1) = obj%rhs_v1(1) * theta0
obj%rhs_v1(2) = obj%rhs_v1(2) * theta0**2

obj%rhs_a1(1) = obj%rhs_a1(2) * theta0**2
obj%rhs_a1(2) = obj%rhs_a1(2) * theta0**3

obj%dis(1) = 1.0_DFP - 1.0_DFP / theta0**3
obj%dis(2) = 1.0_DFP - 1.0_DFP / theta0**2
obj%dis(3) = (1.0_DFP - 1.0_DFP / theta0) * 0.50_DFP
obj%dis(4) = 1.0_DFP / theta0**3

obj%vel(1) = obj%vel(1) / theta0**3
obj%vel(2) = 1.0_DFP - gamma0 / (theta0**2 * beta0)
obj%vel(3) = 1.0_DFP - gamma0 / (2.0_DFP * theta0 * beta0)
obj%vel(4) = obj%vel(4) / theta0**3

obj%acc(1) = obj%acc(1) / theta0**3
obj%acc(2) = obj%acc(2) / theta0**2
obj%acc(3) = 1.0_DFP - 1.0_DFP / (2.0_DFP * theta0 * beta0)
obj%acc(4) = obj%acc(4) / theta0**3

CALL obj%MakeZeros()

END SUBROUTINE obj_Collocation

!----------------------------------------------------------------------------
! Houbolt
!----------------------------------------------------------------------------

SUBROUTINE obj_Houbolt(obj)
CLASS(SDAlgoParam_), INTENT(INOUT) :: obj
CHARACTER(*), PARAMETER :: myName = "Houbolt"
CALL e%RaiseError(modName//'::'//myName//' - '// &
& '[WIP] not implemented yet')
END SUBROUTINE obj_Houbolt

!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
Expand Down Expand Up @@ -229,6 +345,9 @@ SUBROUTINE obj_Display(obj, msg, unitno)
CALL Display(obj%rhs_f2, "rhs_f2: ", unitno=unitno, advance="NO")
CALL Display(obj%rhs_f2_zero, "rhs_f2_zero: ", unitno=unitno)

CALL BlankLines(unitno=unitno)
CALL Display(obj%alpha, "alpha: ", unitno=unitno, advance="NO")

CALL BlankLines(unitno=unitno)
CALL Display(obj%dis, "dis: ", unitno=unitno, advance="NO")
CALL Display(obj%dis_zero, "dis_zero: ", unitno=unitno)
Expand Down Expand Up @@ -274,6 +393,8 @@ SUBROUTINE obj_Deallocate(obj)
obj%vel_zero = .TRUE.
obj%acc_zero = .TRUE.

obj%alpha = 0.0_DFP

END SUBROUTINE obj_Deallocate

!----------------------------------------------------------------------------
Expand All @@ -287,7 +408,7 @@ SUBROUTINE obj_ImportFromToml1(obj, table)
CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml1()"
TYPE(toml_table), POINTER :: node, node2
INTEGER(I4B) :: origin, stat
REAL(DFP) :: beta, gamma
REAL(DFP) :: alpha, beta, gamma, theta, areal
CHARACTER(:), ALLOCATABLE :: str1, algorithm
LOGICAL(LGT) :: problem

Expand Down Expand Up @@ -330,6 +451,48 @@ SUBROUTINE obj_ImportFromToml1(obj, table)
END IF
CALL obj%NewmarkBeta(beta=beta, gamma=gamma)

CASE ("HHTALPHA")

CALL toml_get(table, algorithm, node, origin=origin, requested=.FALSE., &
& stat=stat)

IF (.NOT. ASSOCIATED(node)) THEN
alpha = -0.30_DFP
beta = (1.0_DFP - alpha)**2 * 0.25_DFP
gamma = (1.0_DFP - 2.0_DFP * alpha) * 0.50_DFP
ELSE
CALL toml_get(node, "alpha", alpha, -0.30_DFP, &
& origin=origin, stat=stat)
areal = (1.0_DFP - alpha)**2 * 0.25_DFP
CALL toml_get(node, "beta", beta, areal, &
& origin=origin, stat=stat)
areal = (1.0_DFP - 2.0_DFP * alpha) * 0.50_DFP
CALL toml_get(node, "gamma", gamma, areal, &
& origin=origin, stat=stat)
END IF

CALL obj%HHTAlpha(alpha=alpha, beta=beta, gamma=gamma)

CASE ("COLLOCATION")

CALL toml_get(table, algorithm, node, origin=origin, requested=.FALSE., &
& stat=stat)

IF (.NOT. ASSOCIATED(node)) THEN
beta = 1.0_DFP / 6.0_DFP
gamma = 0.50_DFP
theta = 1.4_DFP
ELSE
CALL toml_get(node, "beta", beta, 1.0_DFP / 6.0_DFP, &
& origin=origin, stat=stat)
CALL toml_get(node, "gamma", gamma, 0.50_DFP, &
& origin=origin, stat=stat)
CALL toml_get(node, "theta", theta, 1.4_DFP, &
& origin=origin, stat=stat)
END IF

CALL obj%Collocation(beta=beta, gamma=gamma, theta=theta)

CASE DEFAULT
node => NULL()
CALL toml_get(table, algorithm, node, origin=origin, requested=.FALSE., &
Expand Down Expand Up @@ -361,6 +524,10 @@ SUBROUTINE obj_ImportFromToml1(obj, table)

IF (.NOT. ASSOCIATED(node2)) THEN
obj%rhs_u1 = 0.0_DFP
obj%rhs_v1 = 0.0_DFP
obj%rhs_a1 = 0.0_DFP
obj%rhs_f1 = 0.0_DFP
obj%rhs_f2 = 0.0_DFP
ELSE

CALL toml_get(node2, "MU", obj%rhs_u1(1), 0.0_DFP, origin=origin, &
Expand Down Expand Up @@ -397,6 +564,9 @@ SUBROUTINE obj_ImportFromToml1(obj, table)
& stat=stat)
END IF

CALL toml_get(node, "alpha", obj%alpha, 0.0_DFP, origin=origin, &
& stat=stat)

CALL toml_get(node, "dis", node2, origin=origin, requested=.FALSE., &
& stat=stat)

Expand Down

0 comments on commit 50e0067

Please sign in to comment.