From aa540cf1326cc99bdf90bd2c3c441288380fbe00 Mon Sep 17 00:00:00 2001 From: shion Date: Tue, 4 Nov 2025 12:17:26 +0900 Subject: [PATCH 1/2] backup --- ...SDAlgorithms.F90 => old_CSDAlgorithms.F90} | 0 src/modules/CMakeLists.txt | 4 + src/modules/CSDAlgorithms/CMakeLists.txt | 21 + .../CSDAlgorithms/src/CSDAlgorithms.F90 | 587 ++++++++++ src/modules/GnuPlot/src/GnuPlot_Class.F90 | 1 + src/modules/SDAlgorithms/src/SDAlgorithms.F90 | 12 +- src/modules/Toml/src/TomlUtility.F90 | 30 +- src/modules/UP2DFEM/CMakeLists.txt | 19 + src/modules/UP2DFEM/src/UP2DFEM_Class.F90 | 1006 +++++++++++++++++ ...tractLinSolver_Class@ImportTomlMethods.F90 | 24 +- .../CSVFile/src/CSVFile_Class@ReadMethods.F90 | 2 +- .../FEDOF/src/FEDOF_Class@GetMethods.F90 | 56 +- .../src/GnuPlot_Class@ConstructorMethods.F90 | 13 + .../GnuPlot/src/GnuPlot_Class@SetMethods.F90 | 1 + .../LinSolver/src/LIS_SOLVE_PRECOND.F90 | 67 ++ .../src/LinSolver_Class@SolveMethods.F90 | 102 +- .../MatrixField_Class@ConstructorMethods.F90 | 36 +- .../Toml/src/TomlUtility@GetMethods.F90 | 64 ++ src/submodules/UP2DFEM/CMakeLists.txt | 29 + .../UP2DFEM_Class@ApplyDirichletBCMethods.F90 | 65 ++ .../src/UP2DFEM_Class@AssembleMethods.F90 | 374 ++++++ .../src/UP2DFEM_Class@ConstructorMethods.F90 | 327 ++++++ .../UP2DFEM/src/UP2DFEM_Class@IOMethods.F90 | 495 ++++++++ .../UP2DFEM/src/UP2DFEM_Class@RunMethods.F90 | 57 + .../UP2DFEM/src/UP2DFEM_Class@SetMethods.F90 | 154 +++ .../src/UP2DFEM_Class@SolveMethods.F90 | 51 + .../src/UP2DFEM_Class@UpdateMethods.F90 | 152 +++ .../src/UP2DFEM_Class@WriteDataMethods.F90 | 91 ++ 28 files changed, 3782 insertions(+), 58 deletions(-) rename src/modules/AbstractKernel/src/{CSDAlgorithms.F90 => old_CSDAlgorithms.F90} (100%) create mode 100644 src/modules/CSDAlgorithms/CMakeLists.txt create mode 100644 src/modules/CSDAlgorithms/src/CSDAlgorithms.F90 create mode 100644 src/modules/UP2DFEM/CMakeLists.txt create mode 100644 src/modules/UP2DFEM/src/UP2DFEM_Class.F90 create mode 100644 src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 create mode 100644 src/submodules/UP2DFEM/CMakeLists.txt create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@ApplyDirichletBCMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@AssembleMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@ConstructorMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@IOMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@RunMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@SetMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@SolveMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@UpdateMethods.F90 create mode 100644 src/submodules/UP2DFEM/src/UP2DFEM_Class@WriteDataMethods.F90 diff --git a/src/modules/AbstractKernel/src/CSDAlgorithms.F90 b/src/modules/AbstractKernel/src/old_CSDAlgorithms.F90 similarity index 100% rename from src/modules/AbstractKernel/src/CSDAlgorithms.F90 rename to src/modules/AbstractKernel/src/old_CSDAlgorithms.F90 diff --git a/src/modules/CMakeLists.txt b/src/modules/CMakeLists.txt index c73a43103..277bd477e 100644 --- a/src/modules/CMakeLists.txt +++ b/src/modules/CMakeLists.txt @@ -408,3 +408,7 @@ include(${CMAKE_CURRENT_LIST_DIR}/Plot/CMakeLists.txt) # semi-discrete algorithms include(${CMAKE_CURRENT_LIST_DIR}/SDAlgorithms/CMakeLists.txt) + +# composite semi-discrete algorithms +include(${CMAKE_CURRENT_LIST_DIR}/CSDAlgorithms/CMakeLists.txt) + diff --git a/src/modules/CSDAlgorithms/CMakeLists.txt b/src/modules/CSDAlgorithms/CMakeLists.txt new file mode 100644 index 000000000..862003a69 --- /dev/null +++ b/src/modules/CSDAlgorithms/CMakeLists.txt @@ -0,0 +1,21 @@ +# This program is a part of EASIFEM library Copyright (C) 2020-2021 Vikas +# Sharma, Ph.D +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see +# + +set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/") +target_sources( + ${PROJECT_NAME} + PRIVATE ${src_path}/CSDAlgorithms.F90) diff --git a/src/modules/CSDAlgorithms/src/CSDAlgorithms.F90 b/src/modules/CSDAlgorithms/src/CSDAlgorithms.F90 new file mode 100644 index 000000000..9d632befc --- /dev/null +++ b/src/modules/CSDAlgorithms/src/CSDAlgorithms.F90 @@ -0,0 +1,587 @@ +! This program is a part of EASIFEM library +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +MODULE CSDAlgorithms +USE ApproxUtility, ONLY: OPERATOR(.APPROXEQ.) + +USE Display_Method, ONLY: Display, & + EqualLine, & + BlankLines, & + tostring + +USE ExceptionHandler_Class, ONLY: e + +USE GlobalData, ONLY: DFP, I4B, LGT, & + CHAR_LF, stdout + +USE InputUtility, ONLY: Input + +USE SDAlgorithms, ONLY: SDAlgoParam_ + +USE StringUtility, ONLY: UpperCase + +USE TomlUtility, ONLY: GetValue + +USE TxtFile_Class, ONLY: TxtFile_ + +USE tomlf, ONLY: toml_table, & + & toml_serialize, & + & toml_get => get_value, & + & toml_len => len, & + & toml_array, & + & toml_stat + +IMPLICIT NONE + +PRIVATE +PUBLIC :: CSDAlgoParam_ +CHARACTER(*), PARAMETER :: modName = "CSDAcousticAlgorithms()" + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-01-19 +! summary: declare of composite type sd algorithm + +TYPE :: CSDAlgoParam_ + TYPE(SDAlgoParam_) :: subAlgoParams(2) + + REAL(DFP) :: splitRatio = 0.5_DFP + ! split ratio for a time interval + ! for example, in bathe method, + ! this value corresponds to gamma + + REAL(DFP) :: rhs_subsol(3, 2) = 0.0_DFP + LOGICAL(LGT) :: rhs_subsol_zero(3, 2) = .TRUE. + + REAL(DFP) :: rhs_subf(2) = 0.0_DFP + LOGICAL(LGT) :: rhs_subf_zero(2) = .TRUE. + + REAL(DFP) :: dis_subsol = 0.0_DFP + LOGICAL(LGT) :: dis_subsol_zero = .TRUE. + + REAL(DFP) :: vel_subsol = 0.0_DFP + LOGICAL(LGT) :: vel_subsol_zero = .TRUE. + + REAL(DFP) :: acc_subsol = 0.0_DFP + LOGICAL(LGT) :: acc_subsol_zero = .TRUE. + + LOGICAL(LGT) :: singleStep = .FALSE. + +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) :: Bathe => obj_Bathe_master + PROCEDURE, PASS(obj) :: ImportFromToml1 => obj_ImportFromToml1 + PROCEDURE, PASS(obj) :: ImportFromToml2 => obj_ImportFromToml2 + GENERIC, PUBLIC :: ImportFromToml => ImportFromToml1, ImportFromToml2 +END TYPE CSDAlgoParam_ + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +CONTAINS + +!---------------------------------------------------------------------------- +! GetBParams_gammaBathe +!---------------------------------------------------------------------------- + +SUBROUTINE GetBParams_gammaBathe(bParams, gamma) + REAL(DFP), INTENT(INOUT) :: bParams(3) + REAL(DFP), OPTIONAL, INTENT(in) :: gamma + + CHARACTER(*), PARAMETER :: myName = "GetBParams_gammaBathe" + REAL(DFP) :: gamma0, areal + + gamma0 = Input(default=0.50_DFP, option=gamma) + + IF (gamma0.approxeq.2.0_DFP) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[CONFIG ERROR] gamma = 2.0 is not allowed') + RETURN + END IF + + areal = 2.0_DFP - gamma0 + bParams(1) = 0.50_DFP / areal + bParams(2) = bParams(1) + bParams(3) = (1.0_DFP - gamma0) / areal + +END SUBROUTINE GetBParams_gammaBathe + +!---------------------------------------------------------------------------- +! GetBParams_betaBathe +!---------------------------------------------------------------------------- + +SUBROUTINE GetBParams_betaBathe(bParams, gamma, beta1, beta2) + REAL(DFP), INTENT(INOUT) :: bParams(3) + REAL(DFP), OPTIONAL, INTENT(in) :: gamma + REAL(DFP), OPTIONAL, INTENT(in) :: beta1 + REAL(DFP), OPTIONAL, INTENT(in) :: beta2 + + CHARACTER(*), PARAMETER :: myName = "GetBParams_betaBathe" + REAL(DFP) :: gamma0, beta1_0, beta2_0 + + gamma0 = Input(default=0.50_DFP, option=gamma) + beta1_0 = Input(default=1.0_DFP / 3.0_DFP, option=beta1) + beta2_0 = Input(default=1.0_DFP - beta1_0, option=beta2) + + bParams(1) = (1.0_DFP - beta1_0) * gamma0 + bParams(2) = (1.0_DFP - beta2_0) * (1.0_DFP - gamma0) + beta1_0 * gamma0 + bParams(3) = (1.0_DFP - gamma0) * beta2_0 + +END SUBROUTINE GetBParams_betaBathe + +!---------------------------------------------------------------------------- +! GetBParams_rhoBathe +!---------------------------------------------------------------------------- + +SUBROUTINE GetBParams_rhoBathe(bParams, gamma, rhoInf) + REAL(DFP), INTENT(INOUT) :: bParams(3) + REAL(DFP), OPTIONAL, INTENT(in) :: gamma + REAL(DFP), OPTIONAL, INTENT(in) :: rhoInf + + CHARACTER(*), PARAMETER :: myName = "GetBParams_rhoBathe" + REAL(DFP) :: gamma0, rhoInf0, areal + + gamma0 = Input(default=0.50_DFP, option=gamma) + rhoInf0 = Input(default=0.0_DFP, option=rhoInf) + + IF (rhoInf0 .LT. 0.0_DFP) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[CONFIG ERROR] asymptotic spectral radius must be positive') + RETURN + END IF + + areal = 2.0_DFP * gamma0 * (rhoInf0 - 1.0_DFP) + 4.0_DFP + bParams(2) = (rhoInf0 + 1.0_DFP) / areal + bParams(1) = 0.50_DFP + (gamma0 - 1.0_DFP) * bParams(2) + bParams(3) = 0.50_DFP - gamma0 * bParams(2) + +END SUBROUTINE GetBParams_rhoBathe + +!---------------------------------------------------------------------------- +! GetBParams_master +!---------------------------------------------------------------------------- + +SUBROUTINE GetBParams_master(bParams, gamma, rhoInf, beta1, beta2) + REAL(DFP), INTENT(INOUT) :: bParams(3) + REAL(DFP), OPTIONAL, INTENT(in) :: gamma + REAL(DFP), OPTIONAL, INTENT(in) :: rhoInf + REAL(DFP), OPTIONAL, INTENT(in) :: beta1 + REAL(DFP), OPTIONAL, INTENT(in) :: beta2 + + CHARACTER(*), PARAMETER :: myName = "GetBParams_master" + LOGICAL(LGT) :: isGamma, isRhoInf, isBeta1, isBeta2, problem + + isGamma = PRESENT(gamma) + isRhoInf = PRESENT(rhoInf) + isBeta1 = PRESENT(beta1) + isBeta2 = PRESENT(beta2) + + problem = isGamma .AND. isRhoInf .AND. isBeta1 .AND. isBeta2 + + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[CONFIG ERROR] too much parameters are passed') + RETURN + END IF + + IF (isRhoInf) THEN + CALL GetBParams_rhoBathe(bParams, gamma=gamma, rhoInf=rhoInf) + RETURN + END IF + + IF (isBeta1 .OR. isBeta2) THEN + CALL GetBParams_rhoBathe(bParams, gamma=gamma, rhoInf=rhoInf) + RETURN + END IF + + IF (isGamma) THEN + CALL GetBParams_gammaBathe(bParams, gamma=gamma) + RETURN + END IF + + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[CONFIG ERROR] no case is found ') + +END SUBROUTINE GetBParams_master + +!---------------------------------------------------------------------------- +! Bathe_master +!---------------------------------------------------------------------------- + +SUBROUTINE obj_Bathe_master(obj, gamma, rhoInf, beta1, beta2) + CLASS(CSDAlgoParam_), INTENT(INOUT) :: obj + REAL(DFP), OPTIONAL, INTENT(IN) :: gamma + REAL(DFP), OPTIONAL, INTENT(IN) :: rhoInf + REAL(DFP), OPTIONAL, INTENT(IN) :: beta1 + REAL(DFP), OPTIONAL, INTENT(IN) :: beta2 + REAL(DFP) :: bParams(3) + + ! internal varibales + CHARACTER(*), PARAMETER :: myName = "obj_Bathe_master" + REAL(DFP) :: gamma0, gamma2, b0, b1, b2 + + CALL obj%DEALLOCATE() + + gamma0 = Input(default=0.5_DFP, option=gamma) + + obj%splitRatio = gamma0 + gamma2 = gamma0 * gamma0 + + !! sub-step 1 + !! equivalent to Trapezoidal rule + obj%subAlgoParams(1)%tanmat(1) = 1.0_DFP + obj%subAlgoParams(1)%tanmat(2) = gamma0 * 0.50_DFP + obj%subAlgoParams(1)%tanmat(3) = gamma2 * 0.25_DFP + + obj%rhs_subf(1) = gamma2 * 0.25_DFP + + obj%subAlgoParams(1)%rhs_u1(1) = 1.0_DFP + obj%subAlgoParams(1)%rhs_u1(2) = 0.50_DFP * gamma0 + + obj%subAlgoParams(1)%rhs_v1(1) = gamma0 + obj%subAlgoParams(1)%rhs_v1(2) = gamma2 * 0.25_DFP + + obj%subAlgoParams(1)%rhs_a1(1) = gamma2 * 0.25_DFP + + !! sub-step 2 + CALL GetBParams_master(bParams, gamma=gamma0, rhoInf=rhoInf, & + & beta1=beta1, beta2=beta2) + b0 = bParams(1) + b1 = bParams(2) + b2 = bParams(3) + + obj%subAlgoParams(2)%tanmat(1) = 1.0_DFP + obj%subAlgoParams(2)%tanmat(2) = b2 + obj%subAlgoParams(2)%tanmat(3) = b2**2 + + obj%subAlgoParams(2)%rhs_f2 = b2**2 + + obj%subAlgoParams(2)%rhs_u1(1) = -2.0_DFP * b1 / gamma0 & + & - 4.0_DFP * b1 * b2 / gamma2 + 1.0_DFP + obj%subAlgoParams(2)%rhs_u1(2) = -2.0_DFP * b1 * b2 / gamma0 + b2 + + obj%rhs_subsol(1, 2) = 2.0_DFP * b1 / gamma0 + & + & 4.0_DFP * b1 * b2 / gamma2 + obj%rhs_subsol(2, 2) = 2.0_DFP * b1 * b2 / gamma0 + + obj%subAlgoParams(2)%rhs_v1(1) = -b0 + b1 + b2 - & + & 4.0_DFP * b1 * b2 / gamma0 + obj%subAlgoParams(2)%rhs_v1(2) = -b2 * (b1 - b0) + + obj%subAlgoParams(2)%rhs_a1(1) = -b2 * (b1 - b0) + + obj%subAlgoParams(2)%acc(1) = 2.0_DFP * b1 / (b2**2 * gamma0) & + & - 1.0_DFP / b2**2 + 4.0_DFP * b1 / (b2 * gamma2) + obj%subAlgoParams(2)%acc(2) = (b1 - b0) / b2**2 - 1.0_DFP / b2 & + & + 4.0_DFP * b1 / (b2 * gamma0) + obj%subAlgoParams(2)%acc(3) = (b1 - b0) / b2 + obj%subAlgoParams(2)%acc(4) = 1.0_DFP / b2**2 + + obj%acc_subsol = -2.0_DFP*b1 / (b2**2*gamma0) - 4.0_DFP*b1 /(b2*gamma2) + + obj%subAlgoParams(2)%vel(1) = 2.0_DFP * b1 / (b2 * gamma0) - 1.0_DFP / b2 + obj%subAlgoParams(2)%vel(2) = (b1 - b0) / b2 + obj%subAlgoParams(2)%vel(3) = 0.0_DFP + obj%subAlgoParams(2)%vel(4) = 1.0_DFP / b2 + + obj%vel_subsol = -2.0_DFP * b1 / (b2 * gamma0) + + obj%subAlgoParams(2)%dis(4) = 1.0_DFP + + obj%dis_subsol = 0.0_DFP + + CALL obj%MakeZeros() + +END SUBROUTINE obj_Bathe_master + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE obj_MakeZeros(obj) + CLASS(CSDAlgoParam_), INTENT(INOUT) :: obj + ! internal varibales + REAL(DFP), PARAMETER :: myzero = 0.0_DFP + INTEGER(I4B) :: ii + + DO ii = 1, SIZE(obj%subAlgoParams) + CALL obj%subAlgoParams(ii)%MakeZeros() + END DO + + obj%rhs_subsol_zero = obj%rhs_subsol.approxeq.myzero + obj%rhs_subf_zero = obj%rhs_subf.approxeq.myzero + obj%dis_subsol_zero = obj%dis_subsol.approxeq.myzero + obj%vel_subsol_zero = obj%vel_subsol.approxeq.myzero + obj%acc_subsol_zero = obj%acc_subsol.approxeq.myzero + +END SUBROUTINE obj_MakeZeros + +!---------------------------------------------------------------------------- +! Display +!---------------------------------------------------------------------------- + +SUBROUTINE obj_Display(obj, msg, unitno) + CLASS(CSDAlgoParam_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: msg + INTEGER(I4B), OPTIONAL, INTENT(IN) :: unitno + INTEGER(I4B) :: ii + + CALL Display(msg, unitno=unitno) + CALL EqualLine(unitno=unitno) + CALL BlankLines(unitno=unitno) + + DO ii = 1, SIZE(obj%subAlgoParams) + CALL obj%subAlgoParams(ii)%Display(msg="sub-step "//tostring(ii), & + & unitno=unitno) + END DO + + CALL BlankLines(unitno=unitno) + CALL Display(obj%splitRatio, "splitRatio: ", unitno=unitno, advance="NO") + + CALL BlankLines(unitno=unitno) + CALL Display(obj%rhs_subsol, "rhs_subsol: ", unitno=unitno, advance="NO") + CALL Display(obj%rhs_subsol_zero, "rhs_subsol_zero: ", unitno=unitno) + + CALL BlankLines(unitno=unitno) + CALL Display(obj%rhs_subf, "rhs_subf: ", unitno=unitno, advance="NO") + CALL Display(obj%rhs_subf_zero, "rhs_subf_zero: ", unitno=unitno) + + CALL BlankLines(unitno=unitno) + CALL Display(obj%dis_subsol, "dis_subsol: ", unitno=unitno, advance="NO") + CALL Display(obj%dis_subsol_zero, "dis_subsol_zero: ", unitno=unitno) + + CALL BlankLines(unitno=unitno) + CALL Display(obj%vel_subsol, "vel_subsol: ", unitno=unitno, advance="NO") + CALL Display(obj%vel_subsol_zero, "vel_subsol_zero: ", unitno=unitno) + + CALL BlankLines(unitno=unitno) + CALL Display(obj%acc_subsol, "acc_subsol_zero: ", unitno=unitno, & + & advance="NO") + CALL Display(obj%acc_subsol_zero, "acc_subsol_zero: ", unitno=unitno) + +END SUBROUTINE obj_Display + +!---------------------------------------------------------------------------- +! Deallocate +!---------------------------------------------------------------------------- + +SUBROUTINE obj_Deallocate(obj) + CLASS(CSDAlgoParam_), INTENT(INOUT) :: obj + INTEGER(I4B) :: ii + + DO ii = 1, SIZE(obj%subAlgoParams) + CALL obj%subAlgoParams(ii)%DEALLOCATE() + END DO + + obj%splitRatio = 0.0_DFP + + obj%rhs_subsol = 0.0_DFP + obj%rhs_subf = 0.0_DFP + + obj%rhs_subsol_zero = .TRUE. + obj%rhs_subf_zero = .TRUE. + + obj%dis_subsol = 0.0_DFP + obj%vel_subsol = 0.0_DFP + obj%acc_subsol = 0.0_DFP + + obj%dis_subsol_zero = .TRUE. + obj%vel_subsol_zero = .TRUE. + obj%acc_subsol_zero = .TRUE. + + obj%singleStep = .FALSE. + +END SUBROUTINE obj_Deallocate + +!---------------------------------------------------------------------------- +! ImportFromToml1 +!---------------------------------------------------------------------------- + +SUBROUTINE obj_ImportFromToml1(obj, table) + CLASS(CSDAlgoParam_), INTENT(INOUT) :: obj + TYPE(toml_table), INTENT(INOUT) :: table + + CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml1()" + TYPE(toml_table), POINTER :: node + INTEGER(I4B) :: origin, stat + REAL(DFP) :: beta, gamma, rhoInf, beta1, beta2 + CHARACTER(:), ALLOCATABLE :: str1, algorithm + LOGICAL(LGT) :: problem + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START]') +#endif + + CALL obj%DEALLOCATE() + + CALL toml_get(table, "algorithm", algorithm, origin=origin, stat=stat) + + problem = (.NOT. ALLOCATED(algorithm)) .OR. (stat .NE. toml_stat%success) + + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: cannot find algorithm field in toml table. ' & + & //'algorithm specifies the name of algorithm.') + RETURN + END IF + + str1 = UpperCase(algorithm) + + SELECT CASE (str1) + CASE ("GAMMABATHE", "GBATHE", "BATHE") + + CALL toml_get(table, algorithm, node, origin=origin, requested=.FALSE., & + & stat=stat) + + IF (.NOT. ASSOCIATED(node)) THEN + gamma = 0.50_DFP + ELSE + CALL toml_get(node, "gamma", gamma, 0.5_DFP, & + & origin=origin, stat=stat) + END IF + CALL obj%Bathe(gamma=gamma) + + CASE ("BETABATHE", "BBATHE") + + CALL toml_get(table, algorithm, node, origin=origin, requested=.FALSE., & + & stat=stat) + + IF (.NOT. ASSOCIATED(node)) THEN + gamma = 0.50_DFP + beta1 = 1.0_DFP / 3.0_DFP + beta2 = 1.0_DFP - beta1 + ELSE + CALL toml_get(node, "gamma", gamma, 0.5_DFP, & + & origin=origin, stat=stat) + CALL toml_get(node, "beta1", beta1, 1.0_DFP / 3.0_DFP, & + & origin=origin, stat=stat) + CALL toml_get(node, "beta2", beta2, 1.0_DFP - beta1, & + & origin=origin, stat=stat) + END IF + CALL obj%Bathe(gamma=gamma, beta1=beta1, beta2=beta2) + + CASE ("RHOBATHE", "RBATHE") + + CALL toml_get(table, algorithm, node, origin=origin, requested=.FALSE., & + & stat=stat) + + IF (.NOT. ASSOCIATED(node)) THEN + gamma = 0.50_DFP + rhoInf = 0.0_DFP + ELSE + CALL toml_get(node, "gamma", gamma, 0.5_DFP, & + & origin=origin, stat=stat) + CALL toml_get(node, "rhoInf", rhoInf, 0.0_DFP, & + & origin=origin, stat=stat) + END IF + CALL obj%Bathe(gamma=gamma, rhoInf=rhoInf) + + CASE ("NEWMARK", "NEWMARKBETA") + + CALL toml_get(table, algorithm, node, origin=origin, & + requested=.FALSE., stat=stat) + + IF (.NOT. ASSOCIATED(node)) THEN + beta = 0.25_DFP; gamma = 0.5_DFP + ELSE + CALL toml_get(node, "beta", beta, 0.25_DFP, & + & origin=origin, stat=stat) + CALL toml_get(node, "gamma", gamma, 0.5_DFP, & + & origin=origin, stat=stat) + END IF + + CALL obj%subAlgoParams(1)%NewmarkBeta(beta=beta, gamma=gamma) + obj%splitRatio = 1.0_DFP + obj%singleStep = .TRUE. + + CASE DEFAULT + + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[CONFIG ERROR] no case is found ') + + END SELECT + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END]') +#endif + +END SUBROUTINE obj_ImportFromToml1 + +!---------------------------------------------------------------------------- +! ImportFromToml2 +!---------------------------------------------------------------------------- + +SUBROUTINE obj_ImportFromToml2(obj, tomlName, afile, & + & filename, printToml) + CLASS(CSDAlgoParam_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: tomlName + TYPE(TxtFile_), OPTIONAL, INTENT(INOUT) :: afile + CHARACTER(*), OPTIONAL, INTENT(IN) :: filename + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: printToml + ! internal variables + CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml2()" + TYPE(toml_table), ALLOCATABLE :: table + TYPE(toml_table), POINTER :: node + INTEGER(I4B) :: origin, stat + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START]') +#endif + + CALL GetValue(table=table, afile=afile, filename=filename) + + node => NULL() + CALL toml_get(table, tomlName, node, origin=origin, requested=.FALSE., & + & stat=stat) + + IF (.NOT. ASSOCIATED(node)) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[CONFIG ERROR] :: following error occured while reading '// & + & 'the toml file :: cannot find ['//tomlName//"] table in config.") + END IF + + CALL obj%ImportFromToml(table=node) + +#ifdef DEBUG_VER + IF (PRESENT(printToml)) THEN + CALL Display(toml_serialize(node), "toml config = "//CHAR_LF, & + & unitNo=stdout) + END IF +#endif + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END]') +#endif + +END SUBROUTINE obj_ImportFromToml2 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +END MODULE CSDAlgorithms diff --git a/src/modules/GnuPlot/src/GnuPlot_Class.F90 b/src/modules/GnuPlot/src/GnuPlot_Class.F90 index 2aae258d8..7d9de97f0 100644 --- a/src/modules/GnuPlot/src/GnuPlot_Class.F90 +++ b/src/modules/GnuPlot/src/GnuPlot_Class.F90 @@ -62,6 +62,7 @@ MODULE GnuPlot_Class TYPE, EXTENDS(AbstractPlot_) :: GnuPlot_ TYPE(TxtFile_) :: pltfile LOGICAL(LGT) :: execute = .TRUE. + LOGICAL(LGT) :: PAUSE = .FALSE. CHARACTER(:), ALLOCATABLE :: commandline TYPE(Label_) :: tpplottitle TYPE(Label_) :: tpxlabel diff --git a/src/modules/SDAlgorithms/src/SDAlgorithms.F90 b/src/modules/SDAlgorithms/src/SDAlgorithms.F90 index d6e5199f1..015919177 100644 --- a/src/modules/SDAlgorithms/src/SDAlgorithms.F90 +++ b/src/modules/SDAlgorithms/src/SDAlgorithms.F90 @@ -142,7 +142,7 @@ SUBROUTINE obj_NewmarkBeta(obj, beta, gamma) ! internal varibales REAL(DFP) :: beta_inv, gamma_inv, beta0, gamma0 - CALL obj%DEALLOCATE() + ! CALL obj%DEALLOCATE() beta0 = Input(default=0.25_DFP, option=beta) gamma0 = Input(default=0.5_DFP, option=gamma) @@ -197,7 +197,7 @@ SUBROUTINE obj_HHTAlpha(obj, alpha, beta, gamma) ! internal varibales REAL(DFP) :: alpha0, beta0, gamma0, areal - CALL obj%DEALLOCATE() + ! CALL obj%DEALLOCATE() alpha0 = Input(default=-0.30_DFP, option=alpha) areal = (1.0_DFP - alpha0)**2 * 0.25_DFP @@ -240,7 +240,7 @@ SUBROUTINE obj_Collocation(obj, beta, gamma, theta) ! internal varibales REAL(DFP) :: theta0, beta0, gamma0, areal - CALL obj%DEALLOCATE() + ! CALL obj%DEALLOCATE() beta0 = Input(default=1.0_DFP / 6.0_DFP, option=beta) gamma0 = Input(default=0.50_DFP, option=gamma) @@ -459,7 +459,7 @@ SUBROUTINE obj_ImportFromToml1(obj, table) CASE ("HHTA") CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & - & stat=stat) + & stat=stat) IF (.NOT. ASSOCIATED(node)) THEN alpha = -0.30_DFP @@ -481,7 +481,7 @@ SUBROUTINE obj_ImportFromToml1(obj, table) CASE ("COLL") CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & - & stat=stat) + & stat=stat) IF (.NOT. ASSOCIATED(node)) THEN beta = 1.0_DFP / 6.0_DFP @@ -501,7 +501,7 @@ SUBROUTINE obj_ImportFromToml1(obj, table) CASE DEFAULT node => NULL() CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & - & stat=stat) + & stat=stat) IF (.NOT. ASSOCIATED(node)) THEN CALL e%RaiseError(modName//'::'//myName//' - '// & diff --git a/src/modules/Toml/src/TomlUtility.F90 b/src/modules/Toml/src/TomlUtility.F90 index 29a1c7c0f..6a4815a33 100644 --- a/src/modules/Toml/src/TomlUtility.F90 +++ b/src/modules/Toml/src/TomlUtility.F90 @@ -50,8 +50,8 @@ MODULE SUBROUTINE GetValue_bool(table, key, VALUE, default_value, & origin, stat, isFound) TYPE(toml_table), INTENT(INOUT) :: table CHARACTER(*), INTENT(IN) :: key - LOGICAL( LGT ), INTENT(INOUT) :: VALUE - LOGICAL( LGT ), INTENT(IN) :: default_value + LOGICAL(LGT), INTENT(INOUT) :: VALUE + LOGICAL(LGT), INTENT(IN) :: default_value INTEGER(I4B), OPTIONAL, INTENT(INOUT) :: origin INTEGER(I4B), OPTIONAL, INTENT(INOUT) :: stat LOGICAL(LGT), OPTIONAL, INTENT(INOUT) :: isFound @@ -86,6 +86,32 @@ MODULE SUBROUTINE GetValue_string(table, key, VALUE, default_value, & END SUBROUTINE GetValue_string END INTERFACE GetValue +!---------------------------------------------------------------------------- +! GetValue@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-08-03 +! summary: GetValue of vector of String + +INTERFACE GetValue + MODULE SUBROUTINE GetValue_string_r1(table, key, VALUE, & + origin, stat, isFound) + TYPE(toml_table), INTENT(INOUT) :: table + !! Toml table + CHARACTER(*), INTENT(IN) :: key + !! key + TYPE(String), ALLOCATABLE, INTENT(INOUT) :: VALUE(:) + !! value in string + INTEGER(I4B), OPTIONAL, INTENT(INOUT) :: origin + !! origin, necessary for debugging + INTEGER(I4B), OPTIONAL, INTENT(INOUT) :: stat + !! To check the status of getting the value + LOGICAL(LGT), OPTIONAL, INTENT(INOUT) :: isFound + !! If key is found then isFound is set to true + END SUBROUTINE GetValue_string_r1 +END INTERFACE GetValue + !---------------------------------------------------------------------------- ! GetValue@Methods !---------------------------------------------------------------------------- diff --git a/src/modules/UP2DFEM/CMakeLists.txt b/src/modules/UP2DFEM/CMakeLists.txt new file mode 100644 index 000000000..038b3b637 --- /dev/null +++ b/src/modules/UP2DFEM/CMakeLists.txt @@ -0,0 +1,19 @@ +# This program is a part of EASIFEM library Copyright (C) 2020-2021 Vikas +# Sharma, Ph.D +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see +# + +set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/") +target_sources(${PROJECT_NAME} PRIVATE ${src_path}/UP2DFEM_Class.F90) diff --git a/src/modules/UP2DFEM/src/UP2DFEM_Class.F90 b/src/modules/UP2DFEM/src/UP2DFEM_Class.F90 new file mode 100644 index 000000000..43f68c2a3 --- /dev/null +++ b/src/modules/UP2DFEM/src/UP2DFEM_Class.F90 @@ -0,0 +1,1006 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +MODULE UP2DFEM_Class +USE GlobalData, ONLY: DFP, I4B, LGT + +USE ParameterList, ONLY: ParameterList_ + +USE BaseType, ONLY: QuadraturePoint_, & + ElemshapeData_, & + poly => TypePolynomialOpt, & + qp => TypeQuadratureOpt, & + ip => TypeInterpolationOpt, & + CSRMatrix_, & + DOF_, & + RealVector_, & + iface_SpaceTimeFunction, & + iface_1DFunction, & + FEVariable_ + +USE tomlf, ONLY: toml_table, & + toml_serialize, & + toml_get => get_value, & + toml_stat, toml_array, & + toml_len => len + +USE FPL, ONLY: ParameterList_ + +USE TxtFile_Class, ONLY: TxtFile_ + +USE ExceptionHandler_Class, ONLY: e + +USE Display_Method, ONLY: Display, ToString + +USE String_Class, ONLY: String + +USE UserFunction_Class, ONLY: UserFunction_ + +USE GnuPlot_Class, ONLY: GnuPlot_ + +USE CSVFile_Class, ONLY: CSVFile_ + +USE SDAlgorithms, ONLY: SDAlgoParam_ + +USE ReallocateUtility, ONLY: Reallocate + +USE FEDomain_Class, ONLY: FEDomain_ +USE FEDOF_Class, ONLY: FEDOF_, FEDOFPointer_ +USE AbstractMesh_Class, ONLY: AbstractMesh_ +USE MatrixField_Class, ONLY: MatrixField_ +USE VectorField_Class, ONLY: VectorField_ +USE NeumannBC_Class, ONLY: NeumannBCPointer_ +USE DirichletBC_Class, ONLY: DirichletBCPointer_ +USE SolidMaterial_Class, ONLY: SolidMaterialPointer_ +USE LinSolver_Class, ONLY: LinSolver_ +USE FEMesh_Class, ONLY: FEMesh_, FEMeshPointer_ + +PRIVATE + +PUBLIC :: UP2DFEM_ + +CHARACTER(*), PARAMETER :: modName = 'UP2DFEM_Class' +CHARACTER(*), PARAMETER :: prefix = "UP2DFEM" +CHARACTER(*), PARAMETER :: default_result_dir = "./results" +CHARACTER(*), PARAMETER :: default_filename = "UP2DFEM" +CHARACTER(*), PARAMETER :: default_baseInterpolationForSpace = "LAGR" +CHARACTER(*), PARAMETER :: default_timeIntegrationScheme = "NEWM" +CHARACTER(*), PARAMETER :: default_baseTypeForSpace = "Monomial" +CHARACTER(*), PARAMETER :: default_ipTypeForSpace = "Equidistance" +CHARACTER(*), PARAMETER :: default_quadTypeForSpace = "GaussLegendre" +INTEGER(I4B), PARAMETER :: MAX_ORDER_SPACE = 10 +INTEGER(I4B), PARAMETER :: default_verbosity = 0 +INTEGER(I4B), PARAMETER :: nsd = 2 + +REAL(DFP), PARAMETER :: one = 1.0_DFP, zero = 0.0_DFP, minus_one = -1.0_DFP, & + half = 0.5_DFP + +#ifdef DEBUG_VER +LOGICAL(LGT), PARAMETER :: debug = .TRUE. +#else +LOGICAL(LGT), PARAMETER :: debug = .FALSE. +#endif + +INTERFACE ElementDataImportFromToml + MODULE PROCEDURE ElementDataImportFromToml_real + MODULE PROCEDURE ElementDataImportFromToml_int +END INTERFACE ElementDataImportFromToml + +!---------------------------------------------------------------------------- +! UP2DFEM_ +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-25 +! summary: Define UP2DFEM_ + +TYPE :: UP2DFEM_ + TYPE(ParameterList_) :: param + TYPE(SDAlgoParam_) :: algoParam + TYPE(FEDomain_) :: dom + TYPE(FEDOF_) :: fedof, fedofNBC + ! TYPE(FEDOFPointer_), ALLOCATABLE :: fedofNBC(:) + + CLASS(AbstractMesh_), POINTER :: cellmesh => NULL() + TYPE(FEMeshPointer_), ALLOCATABLE :: boundaries(:) + + INTEGER(I4B) :: totalSpaceElements = 0 + + TYPE(String) :: meshfilename + TYPE(String) :: filename + + INTEGER(I4B), ALLOCATABLE :: cellcon(:) + + INTEGER(I4B) :: maxNNE = 0, maxCON = 0 + + TYPE(MatrixField_) :: tanmat + TYPE(MatrixField_) :: massMat, stiffMat + + TYPE(VectorField_) :: rhs, sol, rhs0 + TYPE(VectorField_) :: force1, force2, tmp1 + TYPE(VectorField_) :: u0, v0, a0 + + INTEGER(I4B) :: tnbc = 0, tdbc = 0 + + TYPE(NeumannBCPointer_), ALLOCATABLE :: nbc(:) + TYPE(DirichletBCPointer_), ALLOCATABLE :: dbc(:) + INTEGER(I4B), ALLOCATABLE :: dbcIndx(:) + + TYPE(LinSolver_) :: linsolver + + CHARACTER(2) :: baseContinuityForSpace = "H1" + + CHARACTER(4) :: baseInterpolationForSpace = "LAGR" + + CHARACTER(4) :: timeIntegrationScheme = "NEWM" + + INTEGER(I4B) :: verbosity = 0 + + INTEGER(I4B) :: totalTimeSteps = 0 + !! total time steps + + INTEGER(I4B) :: totalVertexDOFSpace = 0 + !! Total number of degree of freedom in space + + INTEGER(I4B) :: totalEdgeDOFSpace = 0 + !! total number of degree of freedom in space + + INTEGER(I4B) :: baseTypeForSpace = poly%monomial + !! basisType for space + !! It is used for LagrangeInterpolation + + INTEGER(I4B) :: ipTypeForSpace = ip%Equidistance + !! interpolation point type for space + !! It is used for LagrangeInpolation + + INTEGER(I4B) :: quadTypeForSpace = qp%GaussLegendre + !! quadrature type for space + + INTEGER(I4B) :: maxSpaceOrder = 0 + !! maximum space order + + INTEGER(I4B) :: currentTimeStep = 1 + !! current time step + + REAL(DFP) :: currentTime = 0.0_DFP + !! current time + + REAL(DFP) :: timeRange(2) = 0.0_DFP + !! Total length of time domain + + INTEGER(I4B) :: spaceOrder = 1 + !! space order of each element + + REAL(DFP), ALLOCATABLE :: timeStepSize(:) + !! size of each time step + !! the size should be totalTimeSteps + + TYPE(SolidMaterialPointer_), ALLOCATABLE :: solidMaterial(:) + + TYPE(FEVariable_) :: density, rayleighAlpha, rayleighBeta, & + cijkl + + TYPE(String) :: result_dir + !! Result directory name + + TYPE(QuadraturePoint_) :: quadForSpace, quadForSpaceBnd + !! Quadrature points in space + + TYPE(ElemshapeData_) :: elemsdForSpace, elemsdForSpaceBnd + !! Element shape data for space + + REAL(DFP), ALLOCATABLE :: ks(:, :), ms(:, :), cs(:, :) + !! space stiffness matrix + + TYPE(UserFunction_), POINTER :: bodyForce => NULL() + !! body force + + TYPE(UserFunction_), POINTER :: initialVel => NULL() + !! initial condition for velocity + + TYPE(UserFunction_), POINTER :: initialDisp => NULL() + !! initial condition for displacement + + TYPE(UserFunction_), POINTER :: initialAcc => NULL() + !! initial condition for acceleration + + LOGICAL(LGT) :: saveData(4) + !! boolean to decide write the data of + !! diaplacement, velocity, acceleration, and all + + INTEGER(I4B) :: outputFreq = 1 + !! output frequency + +CONTAINS + + PROCEDURE, PUBLIC, PASS(obj) :: Initiate => obj_Initiate + !! Initiate + + PROCEDURE, PUBLIC, PASS(obj) :: DEALLOCATE => obj_Deallocate + !! Deallocate data + + PROCEDURE, PASS(obj) :: ImportFromToml1 => obj_ImportFromToml1 + !! Import data from toml file + + PROCEDURE, NON_OVERRIDABLE, PASS(obj) :: ImportFromToml2 => & + obj_ImportFromToml2 + !! Import data from toml file + + GENERIC, PUBLIC :: ImportFromToml => ImportFromToml1, & + ImportFromToml2 + + PROCEDURE, PUBLIC, PASS(obj) :: Display => obj_Display + !! Display data + + PROCEDURE, PUBLIC, PASS(obj) :: InitiateDomains => obj_InitiateDomains + + PROCEDURE, PUBLIC, PASS(obj) :: InitiateFEDOF => obj_InitiateFEDOF + + PROCEDURE, PUBLIC, PASS(obj) :: InitiateFields => obj_InitiateFields + !! Initiate tangent matrix + + PROCEDURE, PUBLIC, PASS(obj) :: Set => obj_Set + !! set the problem + + PROCEDURE, PUBLIC, PASS(obj) :: SetInitialAcceleration => & + obj_SetInitialAcceleration + !! Set initial acceleration + + PROCEDURE, PUBLIC, PASS(obj) :: SetInitialVelocity => & + obj_SetInitialVelocity + !! Set initial velocity + + PROCEDURE, PUBLIC, PASS(obj) :: SetInitialDisplacement => & + obj_SetInitialDisplacement + !! Set initial displacement + + PROCEDURE, PUBLIC, PASS(obj) :: AssembleTanmat => obj_AssembleTanmat + !! Assemble Tangent matrix + + PROCEDURE, PUBLIC, PASS(obj) :: AssembleRHS => obj_AssembleRHS + !! Assemble RHS matrix + + PROCEDURE, PUBLIC, PASS(obj) :: AssembleSurfaceSource => obj_AssembleSurfaceSource + !! Assemble RHS matrix + + PROCEDURE, PUBLIC, PASS(obj) :: ApplyDirichletBC => & + obj_ApplyDirichletBC + !! Apply Left Dirichlet boundary condition + + PROCEDURE, PUBLIC, PASS(obj) :: Solve => obj_Solve + !! Solve + + PROCEDURE, PUBLIC, PASS(obj) :: Update => obj_Update + !! Update + + PROCEDURE, PUBLIC, PASS(obj) :: WriteData => obj_WriteData + !! Update + + PROCEDURE, PUBLIC, PASS(obj) :: Run => obj_Run + !! Debug mode + +END TYPE UP2DFEM_ + +!---------------------------------------------------------------------------- +! - Initiate@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-11 +! summary: Initiate by param + +INTERFACE + MODULE SUBROUTINE obj_Initiate(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Initiate +END INTERFACE + +!---------------------------------------------------------------------------- +! Deallocate@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-11 +! summary: Deallocate the data + +INTERFACE + MODULE SUBROUTINE obj_Deallocate(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Deallocate +END INTERFACE + +!---------------------------------------------------------------------------- +! Display@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-11 +! summary: Display the data + +INTERFACE + MODULE SUBROUTINE obj_Display(obj, msg, unitno) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: msg + INTEGER(I4B), OPTIONAL, INTENT(IN) :: unitno + END SUBROUTINE obj_Display +END INTERFACE + +!---------------------------------------------------------------------------- +! ImportFromToml@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-11 +! summary: Import data from toml config + +INTERFACE + MODULE SUBROUTINE obj_ImportFromToml1(obj, table) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + TYPE(toml_table), INTENT(INOUT) :: table + END SUBROUTINE obj_ImportFromToml1 +END INTERFACE + +!---------------------------------------------------------------------------- +! ImportFromToml@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-11 +! summary: Import data from toml config + +INTERFACE + MODULE SUBROUTINE obj_ImportFromToml2(obj, tomlName, afile, filename, & + printToml) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: tomlName + TYPE(TxtFile_), OPTIONAL, INTENT(INOUT) :: afile + CHARACTER(*), OPTIONAL, INTENT(IN) :: filename + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: printToml + END SUBROUTINE obj_ImportFromToml2 +END INTERFACE + +!---------------------------------------------------------------------------- +! Set@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-02 +! summary: Set the problem + +INTERFACE + MODULE SUBROUTINE obj_Set(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Set +END INTERFACE + +!---------------------------------------------------------------------------- +! GetMs@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-05 +! summary: Get Ms matrix + +INTERFACE + MODULE SUBROUTINE obj_GetMs(obj, ans, nrow, ncol) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + REAL(DFP), INTENT(INOUT) :: ans(:, :) + INTEGER(I4B), INTENT(OUT) :: nrow, ncol + END SUBROUTINE obj_GetMs +END INTERFACE + +!---------------------------------------------------------------------------- +! GetKs@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-05 +! summary: Get Ks matrix + +INTERFACE + MODULE SUBROUTINE obj_GetKs(obj, ans, nrow, ncol) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + REAL(DFP), INTENT(INOUT) :: ans(:, :) + INTEGER(I4B), INTENT(OUT) :: nrow, ncol + END SUBROUTINE obj_GetKs +END INTERFACE + +!---------------------------------------------------------------------------- +! GetBodyForce@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-05 +! summary: Get body force function + +INTERFACE + MODULE SUBROUTINE obj_GetBodyForce(obj, ans, tsize, spaceElemNum, & + time, anscoeff, scale) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + REAL(DFP), INTENT(INOUT) :: ans(:) + INTEGER(I4B), INTENT(OUT) :: tsize + INTEGER(I4B), INTENT(IN) :: spaceElemNum + REAL(DFP), INTENT(IN) :: time + REAL(DFP), INTENT(IN) :: anscoeff + REAL(DFP), INTENT(IN) :: scale + END SUBROUTINE obj_GetBodyForce +END INTERFACE + +!---------------------------------------------------------------------------- +! GetTractionRight@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-05 +! summary: Get traction on right + +INTERFACE + MODULE SUBROUTINE obj_GetTractionRight(obj, ans, tsize, time, & + anscoeff, scale) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + REAL(DFP), INTENT(INOUT) :: ans(:) + INTEGER(I4B), INTENT(OUT) :: tsize + REAL(DFP), INTENT(IN) :: time + REAL(DFP), INTENT(IN) :: anscoeff + REAL(DFP), INTENT(IN) :: scale + END SUBROUTINE obj_GetTractionRight +END INTERFACE + +!---------------------------------------------------------------------------- +! GetTractionLeft@Methods +!---------------------------------------------------------------------------- + +INTERFACE + MODULE SUBROUTINE obj_GetTractionLeft(obj, ans, tsize, time, & + anscoeff, scale) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + REAL(DFP), INTENT(INOUT) :: ans(:) + INTEGER(I4B), INTENT(OUT) :: tsize + REAL(DFP), INTENT(IN) :: time + REAL(DFP), INTENT(IN) :: anscoeff + REAL(DFP), INTENT(IN) :: scale + END SUBROUTINE obj_GetTractionLeft +END INTERFACE + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +INTERFACE + MODULE SUBROUTINE obj_InitiateDomains(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_InitiateDomains +END INTERFACE + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-02-27 +! summary: initiate dof + +INTERFACE + MODULE SUBROUTINE obj_InitiateFEDOF(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_InitiateFEDOF +END INTERFACE + +!---------------------------------------------------------------------------- +! InitiateFields@Methods +!---------------------------------------------------------------------------- + +INTERFACE + MODULE SUBROUTINE obj_InitiateFields(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_InitiateFields +END INTERFACE + +!---------------------------------------------------------------------------- +! AssembleTanmat@Methods +!---------------------------------------------------------------------------- + +INTERFACE + MODULE SUBROUTINE obj_AssembleTanmat(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_AssembleTanmat +END INTERFACE + +!---------------------------------------------------------------------------- +! AssembleRHS@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-12 +! summary: Assemble RHS + +INTERFACE + MODULE SUBROUTINE obj_AssembleRHS(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_AssembleRHS +END INTERFACE + +!---------------------------------------------------------------------------- +! AssembleRHS@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-03-01 +! summary: Assemble RHS + +INTERFACE + MODULE SUBROUTINE obj_AssembleSurfaceSource(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_AssembleSurfaceSource +END INTERFACE + +!---------------------------------------------------------------------------- +! GetCs@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-05 +! summary: Get Cs matrix + +INTERFACE + MODULE SUBROUTINE obj_GetCs(obj, alpha, beta, ans, nrow, ncol) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + REAL(DFP), INTENT(IN) :: alpha, beta + REAL(DFP), INTENT(INOUT) :: ans(:, :) + INTEGER(I4B), INTENT(OUT) :: nrow, ncol + END SUBROUTINE obj_GetCs +END INTERFACE + +!---------------------------------------------------------------------------- +! Debug@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-05 +! summary: Debug mode + +INTERFACE + MODULE SUBROUTINE obj_Run(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Run +END INTERFACE + +!---------------------------------------------------------------------------- +! GetConnectivity@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-08-09 +! summary: Get connectivity matrix + +INTERFACE + MODULE SUBROUTINE obj_GetConnectivity(obj, spaceElemNum, ans, tsize) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + INTEGER(I4B), INTENT(IN) :: spaceElemNum + INTEGER(I4B), INTENT(INOUT) :: ans(:) + INTEGER(I4B), INTENT(OUT) :: tsize + END SUBROUTINE obj_GetConnectivity +END INTERFACE + +!---------------------------------------------------------------------------- +! ApplyDirichletBC@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-09-19 +! summary: Apply Dirichlet boundary condition + +INTERFACE + MODULE SUBROUTINE obj_ApplyDirichletBC(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_ApplyDirichletBC +END INTERFACE + +!---------------------------------------------------------------------------- +! SetInitialAcceleration@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-09-19 +! summary: Set initial acceleration + +INTERFACE + MODULE SUBROUTINE obj_SetInitialAcceleration(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_SetInitialAcceleration +END INTERFACE + +!---------------------------------------------------------------------------- +! SetInitialVelocity@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-09-19 +! summary: Set initial velocity + +INTERFACE + MODULE SUBROUTINE obj_SetInitialVelocity(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_SetInitialVelocity +END INTERFACE + +!---------------------------------------------------------------------------- +! SetInitialDisplacement@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-09-19 +! summary: Set initial displacement + +INTERFACE + MODULE SUBROUTINE obj_SetInitialDisplacement(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_SetInitialDisplacement +END INTERFACE + +!---------------------------------------------------------------------------- +! Solve@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-09-19 +! summary: Set initial displacement + +INTERFACE + MODULE SUBROUTINE obj_Solve(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Solve +END INTERFACE + +!---------------------------------------------------------------------------- +! Update@Methods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-09-19 +! summary: Set initial displacement + +INTERFACE + MODULE SUBROUTINE obj_Update(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Update +END INTERFACE + +!---------------------------------------------------------------------------- +! WriteData@Methods +!---------------------------------------------------------------------------- + +INTERFACE + MODULE SUBROUTINE obj_WriteData(obj) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + END SUBROUTINE obj_WriteData +END INTERFACE + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +CONTAINS + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE ElementDataImportFromToml_real(table, key, val, telem, isfound) + TYPE(toml_table), INTENT(INOUT) :: table + CHARACTER(*), INTENT(IN) :: key + REAL(DFP), ALLOCATABLE, INTENT(INOUT) :: val(:) + INTEGER(I4B), INTENT(IN) :: telem + LOGICAL(LGT), INTENT(INOUT) :: isfound + + CHARACTER(*), PARAMETER :: myName = "ElementDataImportFromToml()" + TYPE(toml_array), POINTER :: array + INTEGER(I4B) :: stat, origin, ii, tsize, ncol, & + iostat + LOGICAL(LGT) :: isOk + REAL(DFP) :: temp + TYPE(String) :: filename, ext + TYPE(TxtFile_) :: atxtfile + TYPE(CSVFile_) :: acsvfile + CHARACTER(512) :: iomsg + INTEGER(I4B), ALLOCATABLE :: intvec1(:), intvec2(:) + REAL(DFP), ALLOCATABLE :: realvec(:) + + IF (debug) CALL Display(myName//key) + + CALL toml_get(table, key, array, origin=origin, stat=stat, & + requested=.FALSE.) + + IF (ASSOCIATED(array)) THEN + tsize = toml_len(array) + CALL Reallocate(val, tsize) + DO ii = 1, tsize + CALL toml_get(array, ii, val(ii)) + END DO + NULLIFY (array) + isfound = .TRUE. + RETURN + END IF + + CALL toml_get(table, key, temp, origin=origin, stat=stat) + + IF (stat .EQ. toml_stat%success) THEN + CALL Reallocate(val, 1) + val(1) = temp + isfound = .TRUE. + RETURN + END IF + + CALL toml_get(table, key, filename%raw, origin=origin, stat=stat) + + IF (stat .EQ. toml_stat%success) THEN + ext = filename%extension() + SELECT CASE (ext%chars()) + CASE (".txt") + CALL atxtfile%Initiate(filename=filename%Chars(), & + action="READ", status="OLD") + CALL atxtfile%OPEN() + CALL atxtfile%READ(val=val, iostat=iostat, iomsg=iomsg) + + isok = iostat .NE. 0 .AND. (.NOT. atxtfile%isEOF()) + IF (isok) THEN + STOP "error occur " + filename = "" + RETURN + END IF + + CALL atxtfile%DEALLOCATE() + filename = "" + isfound = .TRUE. + RETURN + CASE (".csv") + CALL acsvfile%Initiate(filename=filename%Chars(), & + action="READ", status="OLD", & + delimiter=",") + CALL acsvfile%OPEN() + CALL acsvfile%SetHeaderIndx(1) + CALL acsvfile%READ() + ncol = acsvfile%Getncols() + SELECT CASE (ncol) + CASE (1) ! only value + + CALL acsvfile%get(1, val=realvec) ! value + isOk = ALLOCATED(realvec) + CALL AssertError_(isok, myname, & + "value column does not have data") + CALL reallocate(val, 1) + val(1) = realvec(1) + + CASE (2) ! element number and value + + CALL acsvfile%get(1, val=intvec1) ! element + CALL acsvfile%get(2, val=realvec) ! value + + isok = ALLOCATED(intvec1) .AND. ALLOCATED(realvec) + CALL AssertError_(isok, myname, & + "element column or value column "// & + " does not have data") + isok = SIZE(intvec1) .EQ. telem + CALL AssertError_(isok, myname, "the size of elemen column "// & + " should be the same as the number of elements") + CALL reallocate(val, telem) + DO ii = 1, telem + val(intvec1(ii)) = realvec(ii) + END DO + + CASE (3) ! start, end, value + + CALL acsvfile%get(1, val=intvec1) ! start + CALL acsvfile%get(2, val=intvec2) ! end + CALL acsvfile%get(3, val=realvec) ! value + isok = ALLOCATED(intvec1) .AND. ALLOCATED(realvec) & + .AND. ALLOCATED(intvec2) + CALL AssertError_(isok, myname, & + "start column, end column or value column"// & + " does not have data") + isok = MINVAL(intvec1) .EQ. one + CALL AssertError_(isok, myname, "start must have 1 as a component") + isok = MAXVAL(intvec2) .EQ. telem + CALL AssertError_(isok, myname, "end must have the integer"// & + " which is the same as the number of elements") + + CALL reallocate(val, telem) + DO ii = 1, SIZE(intvec1) + val(intvec1(ii):intvec2(ii)) = realvec(ii) + END DO + + CASE default + CALL AssertError_(.FALSE., myname, & + "wrong number of columns in csv file ") + END SELECT + + CALL acsvfile%DEALLOCATE() + filename = "" + isfound = .TRUE. + RETURN + END SELECT + END IF + + isfound = .FALSE. + +END SUBROUTINE ElementDataImportFromToml_real + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE ElementDataImportFromToml_int(table, key, val, telem, isfound) + TYPE(toml_table), INTENT(INOUT) :: table + CHARACTER(*), INTENT(IN) :: key + INTEGER(I4B), ALLOCATABLE, INTENT(INOUT) :: val(:) + INTEGER(I4B), INTENT(IN) :: telem + LOGICAL(LGT), INTENT(INOUT) :: isfound + + CHARACTER(*), PARAMETER :: myName = "ElementDataImportFromToml()" + TYPE(toml_array), POINTER :: array + INTEGER(I4B) :: stat, origin, ii, tsize, ncol, & + iostat, temp + LOGICAL(LGT) :: isOk + TYPE(String) :: filename, ext + TYPE(TxtFile_) :: atxtfile + TYPE(CSVFile_) :: acsvfile + CHARACTER(512) :: iomsg + INTEGER(I4B), ALLOCATABLE :: intvec1(:), intvec2(:), intvec3(:) + + IF (debug) CALL Display(myName//key) + + CALL toml_get(table, key, array, origin=origin, stat=stat, & + requested=.FALSE.) + + IF (ASSOCIATED(array)) THEN + tsize = toml_len(array) + CALL Reallocate(val, tsize) + DO ii = 1, tsize + CALL toml_get(array, ii, val(ii)) + END DO + NULLIFY (array) + isfound = .TRUE. + RETURN + END IF + + CALL toml_get(table, key, temp, origin=origin, stat=stat) + + IF (stat .EQ. toml_stat%success) THEN + CALL Reallocate(val, 1) + val(1) = temp + isfound = .TRUE. + RETURN + END IF + + CALL toml_get(table, key, filename%raw, origin=origin, stat=stat) + + IF (stat .EQ. toml_stat%success) THEN + ext = filename%extension() + SELECT CASE (ext%chars()) + CASE (".txt") + CALL atxtfile%Initiate(filename=filename%Chars(), & + action="READ", status="OLD") + CALL atxtfile%OPEN() + CALL atxtfile%READ(val=val, iostat=iostat, iomsg=iomsg) + + isok = iostat .NE. 0 .AND. (.NOT. atxtfile%isEOF()) + IF (isok) THEN + STOP "error occur " + filename = "" + RETURN + END IF + + CALL atxtfile%DEALLOCATE() + filename = "" + isfound = .TRUE. + RETURN + CASE (".csv") + CALL acsvfile%Initiate(filename=filename%Chars(), & + action="READ", status="OLD", & + delimiter=",") + CALL acsvfile%OPEN() + CALL acsvfile%SetHeaderIndx(1) + CALL acsvfile%READ() + ncol = acsvfile%Getncols() + SELECT CASE (ncol) + CASE (1) ! only value + + CALL acsvfile%get(1, val=intvec3) ! value + isOk = ALLOCATED(intvec3) + CALL AssertError_(isok, myname, & + "value column does not have data") + CALL reallocate(val, 1) + val(1) = intvec3(1) + + CASE (2) ! element number and value + + CALL acsvfile%get(1, val=intvec1) ! element + CALL acsvfile%get(2, val=intvec3) ! value + + isok = ALLOCATED(intvec1) .AND. ALLOCATED(intvec3) + CALL AssertError_(isok, myname, & + "element column or value column "// & + " does not have data") + isok = SIZE(intvec1) .EQ. telem + CALL AssertError_(isok, myname, "the size of elemen column "// & + " should be the same as the number of elements") + CALL reallocate(val, telem) + DO ii = 1, telem + val(intvec1(ii)) = intvec3(ii) + END DO + + CASE (3) ! start, end, value + + CALL acsvfile%get(1, val=intvec1) ! start + CALL acsvfile%get(2, val=intvec2) ! end + CALL acsvfile%get(3, val=intvec3) ! value + isok = ALLOCATED(intvec1) .AND. ALLOCATED(intvec3) & + .AND. ALLOCATED(intvec2) + CALL AssertError_(isok, myname, & + "start column, end column or value column"// & + " does not have data") + isok = MINVAL(intvec1) .EQ. one + CALL AssertError_(isok, myname, "start must have 1 as a component") + isok = MAXVAL(intvec2) .EQ. telem + CALL AssertError_(isok, myname, "end must have the integer"// & + " which is the same as the number of elements") + + CALL reallocate(val, telem) + DO ii = 1, SIZE(intvec1) + val(intvec1(ii):intvec2(ii)) = intvec3(ii) + END DO + + CASE default + CALL AssertError_(.FALSE., myname, & + "wrong number of columns in csv file ") + END SELECT + + CALL acsvfile%DEALLOCATE() + filename = "" + isfound = .TRUE. + RETURN + END SELECT + END IF + + isfound = .FALSE. + +END SUBROUTINE ElementDataImportFromToml_int + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE AssertError_(a, myName, msg) + LOGICAL(LGT), INTENT(IN) :: a + CHARACTER(*), INTENT(IN) :: myName + CHARACTER(*), INTENT(IN) :: msg + + IF (.NOT. a) THEN + CALL e%RaiseError(modName//'::'//myName//" - "// & + '[INTERNAL ERROR] :: '//msg) + RETURN + END IF + +END SUBROUTINE AssertError_ + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +END MODULE UP2DFEM_Class diff --git a/src/submodules/AbstractLinSolver/src/AbstractLinSolver_Class@ImportTomlMethods.F90 b/src/submodules/AbstractLinSolver/src/AbstractLinSolver_Class@ImportTomlMethods.F90 index b53c29d8e..d50a654a2 100644 --- a/src/submodules/AbstractLinSolver/src/AbstractLinSolver_Class@ImportTomlMethods.F90 +++ b/src/submodules/AbstractLinSolver/src/AbstractLinSolver_Class@ImportTomlMethods.F90 @@ -245,7 +245,8 @@ END SUBROUTINE iluc_import_from_toml CHARACTER(:), ALLOCATABLE :: prefix TYPE(String) :: engine, solverName, preconditionOption, & - p_name, convergenceIn, convergenceType, scale, p_hybrid_i + p_name, convergenceIn, convergenceType, scale, p_hybrid_i, & + astr INTEGER(I4B) :: maxIter, krylovSubspaceSize, bicgstab_ell, & p_ilu_lfil, p_ilu_mbloc, p_ilu_fill, p_hybrid_maxiter, & @@ -387,25 +388,26 @@ END SUBROUTINE iluc_import_from_toml CALL toml_get(node, LowerCase(p_name%chars()), child, origin=origin, & stat=stat, requested=.FALSE.) IF (ASSOCIATED(child)) THEN - SELECT CASE (p_name%chars()) - CASE ("none", "NONE") + astr = p_name%upper() + SELECT CASE (astr%chars()) + CASE ("NONE") ! do nothing - CASE ("ilu", "ILU") + CASE ("ILU", "ILU0", "ILUT", "ILUTP", "ILUK", "ILUD", "ILUDP") CALL ilu_import_from_toml(param=param, prefix=prefix, table=child) - CASE ("hybrid", "HYBRID") + CASE ("HYBRID") CALL hybrid_import_from_toml(obj=obj, param=param, prefix=prefix, & table=child) - CASE ("is", "IS") + CASE ("IS") CALL is_import_from_toml(param=param, prefix=prefix, table=child) - CASE ("adds", "ADDS") + CASE ("ADDS") CALL adds_import_from_toml(param=param, prefix=prefix, table=child) - CASE ("ssor", "SSOR") + CASE ("SSOR") CALL ssor_import_from_toml(param=param, prefix=prefix, table=child) - CASE ("sainv", "SAINV") + CASE ("SAINV") CALL sainv_import_from_toml(param=param, prefix=prefix, table=child) - CASE ("saamg", "SAAMG") + CASE ("SAAMG") CALL saamg_import_from_toml(param=param, prefix=prefix, table=child) - CASE ("iluc", "ILUC") + CASE ("ILUC") CALL iluc_import_from_toml(param=param, prefix=prefix, table=child) END SELECT child => NULL() diff --git a/src/submodules/CSVFile/src/CSVFile_Class@ReadMethods.F90 b/src/submodules/CSVFile/src/CSVFile_Class@ReadMethods.F90 index 3ed632263..12415ca56 100644 --- a/src/submodules/CSVFile/src/CSVFile_Class@ReadMethods.F90 +++ b/src/submodules/CSVFile/src/CSVFile_Class@ReadMethods.F90 @@ -51,7 +51,7 @@ isSkipRows = .FALSE. END IF -nrows = trecords - skippedRows - obj%headerIndx +nrows = trecords - skippedRows - MIN(1_I4B, obj%headerIndx) obj%nrows = nrows trecords = obj%getTotalRecords() diff --git a/src/submodules/FEDOF/src/FEDOF_Class@GetMethods.F90 b/src/submodules/FEDOF/src/FEDOF_Class@GetMethods.F90 index 7363dc8c2..04cc55e1c 100644 --- a/src/submodules/FEDOF/src/FEDOF_Class@GetMethods.F90 +++ b/src/submodules/FEDOF/src/FEDOF_Class@GetMethods.F90 @@ -29,6 +29,8 @@ USE Display_Method, ONLY: ToString +USE StringUtility, ONLY: UpperCase + #ifdef DEBUG_VER USE Display_Method, ONLY: Display #endif @@ -223,20 +225,27 @@ INTEGER(I4B) :: ent(4) INTEGER(I4B) :: ii, jj, kk, a, b INTEGER(I4B) :: temp(PARAM_MAX_CONNECTIVITY_SIZE) +LOGICAL(LGT) :: isok +CHARACTER(1) :: opt0 ent = obj%mesh%GetTotalEntities(globalElement=globalElement, islocal=islocal) CALL obj%mesh%GetConnectivity_(globalElement=globalElement, islocal=islocal, & - opt=opt, tsize=jj, ans=temp) + opt="A", tsize=jj, ans=temp) + +opt0 = Uppercase(opt(1:1)) ! points a = 1; b = ent(1) jj = 1 -DO ii = a, b - CALL obj%GetVertexDOF(globalNode=temp(ii), ans=ans(jj:), tsize=kk, & - islocal=.FALSE.) - jj = jj + kk -END DO +isok = opt0 .EQ. "V" .OR. opt0 .EQ. "A" +IF (isok) THEN + DO ii = a, b + CALL obj%GetVertexDOF(globalNode=temp(ii), ans=ans(jj:), tsize=kk, & + islocal=.FALSE.) + jj = jj + kk + END DO +END IF IF (obj%isLagrange) THEN tsize = jj - 1 @@ -245,25 +254,34 @@ ! edges a = b + 1; b = b + ent(2) -DO ii = a, b - CALL obj%GetEdgeDOF(globalEdge=temp(ii), ans=ans(jj:), tsize=kk) - jj = jj + kk -END DO +isok = opt0 .EQ. "E" .OR. opt0 .EQ. "A" +IF (isok) THEN + DO ii = a, b + CALL obj%GetEdgeDOF(globalEdge=temp(ii), ans=ans(jj:), tsize=kk) + jj = jj + kk + END DO +END IF ! faces a = b + 1; b = b + ent(3) -DO ii = a, b - CALL obj%GetFaceDOF(globalFace=temp(ii), ans=ans(jj:), tsize=kk) - jj = jj + kk -END DO +isok = opt0 .EQ. "F" .OR. opt0 .EQ. "A" +IF (isok) THEN + DO ii = a, b + CALL obj%GetFaceDOF(globalFace=temp(ii), ans=ans(jj:), tsize=kk) + jj = jj + kk + END DO +END IF ! cell a = b + 1; b = b + ent(4) -DO ii = a, b - CALL obj%GetCellDOF(globalCell=temp(ii), ans=ans(jj:), tsize=kk, & - islocal=.FALSE.) - jj = jj + kk -END DO +isok = opt0 .EQ. "C" .OR. opt0 .EQ. "A" +IF (isok) THEN + DO ii = a, b + CALL obj%GetCellDOF(globalCell=temp(ii), ans=ans(jj:), tsize=kk, & + islocal=.FALSE.) + jj = jj + kk + END DO +END IF tsize = jj - 1 diff --git a/src/submodules/GnuPlot/src/GnuPlot_Class@ConstructorMethods.F90 b/src/submodules/GnuPlot/src/GnuPlot_Class@ConstructorMethods.F90 index 06c76e4ce..02a3231b4 100644 --- a/src/submodules/GnuPlot/src/GnuPlot_Class@ConstructorMethods.F90 +++ b/src/submodules/GnuPlot/src/GnuPlot_Class@ConstructorMethods.F90 @@ -68,6 +68,19 @@ obj%plotEngine = PLOT_ENGINE_PLPLOT +IF (obj%PAUSE) THEN + CALL obj%pltfile%WRITE("pause mouse keypress") + CALL obj%pltfile%WRITE("if (exists('MOUSE_CHAR') && MOUSE_CHAR eq 'c' || MOUSE_KEY == -1) {") + CALL obj%pltfile%WRITE(" print 'Exiting...'") + CALL obj%pltfile%WRITE(" quit") + CALL obj%pltfile%WRITE(" } else {") + call obj%pltfile%WRITE(" print 'Press c or q to quit, any other key to refresh'") + CALL obj%pltfile%WRITE(" replot") + CALL obj%pltfile%WRITE(" load '"//obj%txtfilename//"'") + CALL obj%pltfile%WRITE(" }") + obj%commandline = "gnuplot " +END IF + IF (obj%pltfile%IsOpen()) THEN CALL obj%pltfile%DEALLOCATE() obj%hasfileopen = .FALSE. diff --git a/src/submodules/GnuPlot/src/GnuPlot_Class@SetMethods.F90 b/src/submodules/GnuPlot/src/GnuPlot_Class@SetMethods.F90 index aedd61a99..4c5b52b3c 100644 --- a/src/submodules/GnuPlot/src/GnuPlot_Class@SetMethods.F90 +++ b/src/submodules/GnuPlot/src/GnuPlot_Class@SetMethods.F90 @@ -316,6 +316,7 @@ obj%commandline = defaultCommandLine obj%execute = .TRUE. +obj%PAUSE = .FALSE. END PROCEDURE reset_to_defaults diff --git a/src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 b/src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 new file mode 100644 index 000000000..c80b5e7d6 --- /dev/null +++ b/src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 @@ -0,0 +1,67 @@ +SUBROUTINE _SUBROUTINE_NAME(obj, sol, rhs, precond) + CLASS(LinSolver_), TARGET, INTENT(INOUT) :: obj + TYPE(CSRMatrix_), INTENT(IN) :: precond + REAL(DFP), INTENT(INOUT) :: sol(:) + REAL(DFP), INTENT(INOUT) :: rhs(:) + + ! Internal variables + CHARACTER(*), PARAMETER :: myName = _MY_NAME + INTEGER(I4B) :: n + REAL(DFP), ALLOCATABLE :: diag(:) + CLASS(AbstractMatrixField_), POINTER :: amat + LOGICAL(LGT) :: isok + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + + obj%IPAR(1) = 0 + obj%FPAR(11) = 0.0_DFP + CALL obj%GetParam(globalNumRow=n, amat=amat) + obj%IPAR(7) = 1 + + DO + + CALL _LIS_NAME(n, rhs, sol, obj%IPAR, obj%FPAR, obj%W) + + IF (obj%IPAR(1) .GT. 0) THEN + + CALL PERFORM_TASK_PRECOND(amat, y=obj%W(obj%IPAR(9):obj%IPAR(9) + n - 1), & + x=obj%W(obj%IPAR(8):obj%IPAR(8) + n - 1), & + precond=precond, & + ierr=obj%IPAR(1)) + + ELSE IF (obj%IPAR(1) .LT. 0) THEN + + CALL CHECKERROR(IPAR=obj%IPAR, FPAR=obj%FPAR, myName=myName) + EXIT + + ELSE IF (obj%IPAR(1) .EQ. 0) THEN + + CALL obj%SetParam(ierr=obj%ipar(1), iter=obj%ipar(7)) + CALL DisplayConvergence(myName, obj%ipar(7), obj%FPAR) + EXIT + + END IF + + END DO + + ! Initial residual/error norm + + CALL obj%SetParam(error0=obj%fpar(3), tol=obj%fpar(4), & + error=obj%fpar(6), normRes=obj%fpar(5)) + +END SUBROUTINE _SUBROUTINE_NAME + +#ifdef _SUBROUTINE_NAME +#undef _SUBROUTINE_NAME +#endif + +#ifdef _LIS_NAME +#undef _LIS_NAME +#endif + +#ifdef _MY_NAME +#undef _MY_NAME +#endif diff --git a/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 b/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 index fcc08ea88..ea4cb3d43 100644 --- a/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 +++ b/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 @@ -21,7 +21,7 @@ USE BaseType, ONLY: TypeSolverNameOpt -USE CSRMatrix_Method, ONLY: LinSolve +USE CSRMatrix_Method, ONLY: LinSolve, MatVec USE SuperLU_Types, ONLY: yes_no_t @@ -82,6 +82,59 @@ SUBROUTINE PERFORM_TASK(amat, y, x, ierr) END SUBROUTINE PERFORM_TASK +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE PERFORM_TASK_PRECOND(amat, y, x, precond, ierr) + + ! intent of dummy variables + CLASS(AbstractMatrixField_), INTENT(INOUT) :: amat + TYPE(CSRMatrix_), INTENT(IN) :: precond + REAL(DFP), INTENT(INOUT) :: y(:) + REAL(DFP), INTENT(IN) :: x(:) + INTEGER(I4B), INTENT(IN) :: ierr + +#ifdef DEBUG_VER + CHARACTER(*), PARAMETER :: myName = "PERFORM_TASK()" +#endif + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + + SELECT CASE (ierr) + + CASE (1) + + CALL amat%Matvec(y=y, x=x, isTranspose=.FALSE.) + + CASE (2) + + CALL amat%Matvec(y=y, x=x, isTranspose=.TRUE.) + + CASE (3, 5) + + ! LEFT/RIGHT PRECONDITIONER SOLVER + ! The preconditioners are given precond + CALL Matvec(obj=precond, x=x, y=y, isTranspose=.FALSE.) + + CASE (4, 6) + + ! LEFT/RIGHT PRECONDITIONER SOLVER + ! The preconditioners are given precond + CALL Matvec(obj=precond, x=x, y=y, isTranspose=.TRUE.) + + END SELECT + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif + +END SUBROUTINE PERFORM_TASK_PRECOND + !---------------------------------------------------------------------------- ! CHECKERR !---------------------------------------------------------------------------- @@ -332,6 +385,43 @@ END SUBROUTINE DisplayConvergence END PROCEDURE obj_Solve +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Solve_precond +CHARACTER(*), PARAMETER :: myName = "obj_Solve()" +REAL(DFP), POINTER :: rhsvar(:), solvar(:) +INTEGER(I4B) :: info +INTEGER(I4B) :: solverName +LOGICAL(LGT) :: isok +CLASS(AbstractMatrixField_), POINTER :: amat + +CALL obj%GetParam(isInitiated=isok, solverName=solverName, amat=amat) + +CALL AssertError1(isok, myname, 'Linear solver is not initiated!') + +isok = ASSOCIATED(amat) +CALL AssertError1(isok, myname, 'Amat is not associated') + +SELECT CASE (solverName) + +CASE (TypeSolverNameOpt%GMRES) + + rhsvar => rhs%GetPointer() + solvar => sol%GetPointer() + CALL LS_SOLVE_GMRES_precond(obj=obj, sol=solvar, rhs=rhsvar, & + precond=precond) + NULLIFY (rhsvar, solvar) + +CASE DEFAULT + + CALL AssertError1(.FALSE., myName, 'No case found for linear solver') + RETURN +END SELECT + +END PROCEDURE obj_Solve_precond + !---------------------------------------------------------------------------- ! LS_SOLVE_CG !---------------------------------------------------------------------------- @@ -405,6 +495,16 @@ END SUBROUTINE DisplayConvergence #include "./LIS_SOLVE.F90" +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#define _SUBROUTINE_NAME LS_SOLVE_GMRES_PRECOND +#define _LIS_NAME GMRES +#define _MY_NAME "LS_SOLVE_GMRES_PRECOND" + +#include "./LIS_SOLVE_PRECOND.F90" + !---------------------------------------------------------------------------- ! LS_SOLVE_FGMRES !---------------------------------------------------------------------------- diff --git a/src/submodules/MatrixField/src/MatrixField_Class@ConstructorMethods.F90 b/src/submodules/MatrixField/src/MatrixField_Class@ConstructorMethods.F90 index e30bb3a2b..7026c880b 100644 --- a/src/submodules/MatrixField/src/MatrixField_Class@ConstructorMethods.F90 +++ b/src/submodules/MatrixField/src/MatrixField_Class@ConstructorMethods.F90 @@ -122,13 +122,13 @@ RETURN END IF -CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="name", VALUE=name) -CALL Set(obj=sublist, datatype="Char", prefix=prefix, key="engine", VALUE=engine) -CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="comm", & +CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="name", VALUE=name) +CALL Set(obj=param, datatype="Char", prefix=prefix, key="engine", VALUE=engine) +CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="comm", & VALUE=INPUT(option=comm, default=0_I4B)) -CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="local_n", & +CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="local_n", & VALUE=INPUT(option=local_n, default=0_I4B)) -CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="global_n", & +CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="global_n", & VALUE=INPUT(option=global_n, default=0_I4B)) SELECT CASE (name) @@ -143,10 +143,10 @@ RETURN END IF - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="droptol", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="droptol", & VALUE=droptol) - CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="lfil", & + CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="lfil", & VALUE=lfil) CASE (PRECOND_ILUTP) @@ -163,16 +163,16 @@ RETURN END IF - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="droptol", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="droptol", & VALUE=droptol) - CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="lfil", & + CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="lfil", & VALUE=lfil) - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="permtol", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="permtol", & VALUE=permtol) - CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="mbloc", & + CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="mbloc", & VALUE=mbloc) CASE (PRECOND_ILUD) @@ -186,9 +186,9 @@ RETURN END IF - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="droptol", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="droptol", & VALUE=droptol) - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="alpha", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="alpha", & VALUE=alpha) CASE (PRECOND_ILUDP) @@ -205,13 +205,13 @@ RETURN END IF - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="droptol", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="droptol", & VALUE=droptol) - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="alpha", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="alpha", & VALUE=alpha) - CALL Set(obj=sublist, datatype=1.0_DFP, prefix=prefix, key="permtol", & + CALL Set(obj=param, datatype=1.0_DFP, prefix=prefix, key="permtol", & VALUE=permtol) - CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="mbloc", & + CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="mbloc", & VALUE=mbloc) CASE (PRECOND_ILUK) @@ -223,7 +223,7 @@ RETURN END IF - CALL Set(obj=sublist, datatype=1_I4B, prefix=prefix, key="lfil", & + CALL Set(obj=param, datatype=1_I4B, prefix=prefix, key="lfil", & VALUE=lfil) CASE DEFAULT diff --git a/src/submodules/Toml/src/TomlUtility@GetMethods.F90 b/src/submodules/Toml/src/TomlUtility@GetMethods.F90 index 6f4cb3732..4f5b47369 100644 --- a/src/submodules/Toml/src/TomlUtility@GetMethods.F90 +++ b/src/submodules/Toml/src/TomlUtility@GetMethods.F90 @@ -34,6 +34,8 @@ USE CSVFile_Class, ONLY: CSVFile_ +USE String_Class, ONLY: reallocate + IMPLICIT NONE CONTAINS @@ -60,6 +62,68 @@ IF (PRESENT(isFound)) isFound = .FALSE. END PROCEDURE GetValue_string +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE GetValue_string_r1 +CHARACTER(:), ALLOCATABLE :: temp_char +LOGICAL(LGT) :: isok + +TYPE(toml_array), POINTER :: array +INTEGER(I4B) :: tsize, ii, stat0 +LOGICAL(LGT) :: isFound0, isok + +isFound0 = .FALSE. + +!---------------------------------------------------------------------------- +! READ from TOML array +! try to read from the toml array +! the data is given in toml file itself as toml array +!---------------------------------------------------------------------------- +array => NULL() +CALL toml_get(table, key, array, origin=origin, stat=stat0, & + requested=.FALSE.) + +isok = ASSOCIATED(array) + +IF (isok) THEN + tsize = toml_len(array) + CALL Reallocate(VALUE, tsize) + isFound0 = .TRUE. + DO ii = 1, tsize + CALL toml_get(array, ii, temp_char) + VALUE(ii) = temp_char + END DO + + IF (PRESENT(stat)) stat = stat0 + IF (PRESENT(isFound)) isFound = isFound0 + NULLIFY (array) + RETURN +END IF + +!---------------------------------------------------------------------------- +! READ a single value from toml +! In this case length of the vector is 1, this value is given in toml file +!---------------------------------------------------------------------------- + +CALL toml_get(table, key, temp_char, origin=origin, stat=stat0) + +IF (stat0 .EQ. toml_stat%success) THEN + CALL Reallocate(VALUE, 1) + VALUE(1) = temp_char + isFound0 = .TRUE. + IF (PRESENT(isFound)) isFound = isFound0 + IF (PRESENT(stat)) stat = stat0 + RETURN +END IF + +isFound0 = .FALSE. +IF (PRESENT(isFound)) isFound = isFound0 +IF (PRESENT(stat)) stat = stat0 + +END PROCEDURE GetValue_string_r1 + !---------------------------------------------------------------------------- ! Get !---------------------------------------------------------------------------- diff --git a/src/submodules/UP2DFEM/CMakeLists.txt b/src/submodules/UP2DFEM/CMakeLists.txt new file mode 100644 index 000000000..bfe79b834 --- /dev/null +++ b/src/submodules/UP2DFEM/CMakeLists.txt @@ -0,0 +1,29 @@ +# This program is a part of EASIFEM library Copyright (C) 2020-2021 Vikas +# Sharma, Ph.D +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see +# + +set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/") +target_sources( + ${PROJECT_NAME} + PRIVATE ${src_path}/UP2DFEM_Class@ConstructorMethods.F90 + ${src_path}/UP2DFEM_Class@IOMethods.F90 + ${src_path}/UP2DFEM_Class@SetMethods.F90 + ${src_path}/UP2DFEM_Class@AssembleMethods.F90 + ${src_path}/UP2DFEM_Class@ApplyDirichletBCMethods.F90 + ${src_path}/UP2DFEM_Class@UpdateMethods.F90 + ${src_path}/UP2DFEM_Class@SolveMethods.F90 + ${src_path}/UP2DFEM_Class@WriteDataMethods.F90 + ${src_path}/UP2DFEM_Class@RunMethods.F90) diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@ApplyDirichletBCMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@ApplyDirichletBCMethods.F90 new file mode 100644 index 000000000..34d05ca1c --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@ApplyDirichletBCMethods.F90 @@ -0,0 +1,65 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) ApplyDirichletBCMethods + +USE Display_Method + +IMPLICIT NONE + +REAL(DFP), PARAMETER :: one = 1.0_DFP, zero = 0.0_DFP, minus_one = -1.0_DFP, & + half = 0.5_DFP + +CONTAINS + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplyDirichletBC +CHARACTER(*), PARAMETER :: myName = "obj_ApplyDirichletBC()" +REAL(DFP) :: times(1) +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +isok = ALLOCATED(obj%dbc) + +IF (isok) THEN + times = obj%currentTime + & + obj%timeStepSize(obj%currentTimeStep) + CALL obj%sol%Set(VALUE=zero) + CALL obj%sol%ApplyDirichletBC(dbc=obj%dbc, times=times) + CALL obj%tanmat%ApplyDBCtoRHS(x=obj%sol, y=obj%rhs, & + scale=minus_one, addContribution=.TRUE.) + CALL obj%rhs%ApplyDirichletBC(dbc=obj%dbc, times=times) +END IF + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_ApplyDirichletBC + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE ApplyDirichletBCMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@AssembleMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@AssembleMethods.F90 new file mode 100644 index 000000000..1632005ba --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@AssembleMethods.F90 @@ -0,0 +1,374 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) AssembleMethods + +USE BaseType, ONLY: TypeFEVariableSpace, & + TypeFEVariableVector, & + TypeFEVariableScalar + +USE ProductUtility, ONLY: OuterProd_, OTimesTilda + +USE GlobalData, ONLY: DOF_FMT, & + NODES_FMT, & + NONE + +USE ReallocateUtility, ONLY: Reallocate + +USE MassMatrix_Method +USE StiffnessMatrix_Method +USE BaseType, ONLY: FEVariableScalar_ +USE NeumannBC_Class, ONLY: NeumannBC_ +USE FEVariable_Method, ONLY: NodalVariable, & + FEVar_Display => Display, & + DEALLOCATE +USE ForceVector_Method, ONLY: ForceVector +USE Display_Method +USE ElemshapeData_InterpolMethods + +IMPLICIT NONE + +CONTAINS + +!---------------------------------------------------------------------------- +! AssembleTanmat +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_AssembleTanmat +CHARACTER(*), PARAMETER :: myName = "obj_AssembleTanmat()" +INTEGER(I4B) :: nrow, ncol, tsize, iel +REAL(DFP) :: dt, dts, scale +REAL(DFP) :: xij(3, obj%maxNNE) +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +dt = obj%timeStepSize(obj%currentTimeStep) +dts = dt * dt + +CALL obj%massMat%set(VALUE=zero) +! CALL obj%dampMat%set(VALUE=zero) +CALL obj%stiffMat%set(VALUE=zero) + +DO iel = 1, obj%totalSpaceElements + + CALL obj%fedof%GetQuadraturePoints1(quad=obj%quadForSpace, & + globalElement=iel, & + quadratureType=obj%quadTypeForSpace, & + order=2 * obj%spaceOrder, & + islocal=.TRUE.) + + CALL obj%fedof%GetLocalElemShapeData(globalElement=iel, & + elemsd=obj%elemsdForSpace, & + quad=obj%quadForSpace, & + islocal=.TRUE.) + + CALL obj%cellmesh%GetNodeCoord(nodecoord=xij, nrow=nrow, ncol=ncol, & + globalElement=iel, islocal=.TRUE.) + + CALL obj%fedof%GetGlobalElemShapeData(globalElement=iel, & + elemsd=obj%elemsdForSpace, & + xij=xij, islocal=.TRUE.) + + CALL obj%fedof%GetConnectivity_(globalElement=iel, islocal=.TRUE., & + ans=obj%cellcon, tsize=tsize, opt="A") + + obj%ms = zero + CALL MassMatrix_(test=obj%elemsdForSpace, trial=obj%elemsdForSpace, & + rho=obj%density, rhorank=TypeFEVariableScalar, & + ans=obj%ms, nrow=nrow, ncol=ncol, opt=nsd) + + CALL obj%massMat%Set(globalNode=obj%cellcon(1:tsize), islocal=.TRUE., & + VALUE=obj%ms(1:nrow, 1:ncol), storageFMT=DOF_FMT, & + scale=one, addContribution=.TRUE.) + + obj%ks = zero + CALL StiffnessMatrix_(test=obj%elemsdForSpace, trial=obj%elemsdForSpace, & + Cijkl=obj%Cijkl, ans=obj%ks, nrow=nrow, ncol=ncol) + + CALL obj%stiffMat%Set(globalNode=obj%cellcon(1:tsize), islocal=.TRUE., & + VALUE=obj%ks(1:nrow, 1:ncol), storageFMT=DOF_FMT, & + scale=one, addContribution=.TRUE.) + +END DO + +CALL obj%tanmat%set(VALUE=zero) + +scale = obj%algoParam%tanmat(1) +CALL obj%tanmat%Set(VALUE=obj%massMat, scale=scale, & + addContribution=.TRUE.) + +scale = obj%algoParam%tanmat(3) * dts +CALL obj%tanmat%Set(VALUE=obj%stiffMat, scale=scale, & + addContribution=.TRUE.) + +isok = ALLOCATED(obj%dbcIndx) +CALL obj%tanmat%ApplyDBC(dbcPtrs=obj%dbcIndx) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_AssembleTanmat + +!---------------------------------------------------------------------------- +! AssembleRHS +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_AssembleRHS + +CHARACTER(*), PARAMETER :: myName = "obj_AssembleRHS()" +REAL(DFP) :: scale, dt, dt2 +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') + +dt = obj%timeStepSize(obj%currentTimeStep) +dt2 = dt * dt + +CALL obj%rhs%Set(VALUE=zero) + +isok = .NOT. obj%algoParam%rhs_f1_zero +IF (isok) THEN + scale = obj%algoParam%rhs_f1 * dt2 + CALL obj%rhs%AXPY(x=obj%force1, scale=scale) +END IF + +isok = .NOT. obj%algoParam%rhs_f2_zero +IF (isok) THEN + scale = obj%algoParam%rhs_f2 * dt2 + CALL obj%rhs%AXPY(x=obj%force2, scale=scale) +END IF + +!---------------------------------------------------------------------------- +! make M(a*U+b*V*dt+c*A*dt*dt) +!---------------------------------------------------------------------------- + +CALL obj%tmp1%Set(VALUE=zero) +scale = obj%algoParam%rhs_u1(1) +isok = .NOT. obj%algoParam%rhs_u1_zero(1) +IF (isok) & + & CALL obj%tmp1%AXPY(x=obj%u0, scale=scale) + +scale = obj%algoParam%rhs_v1(1) * dt +isok = .NOT. obj%algoParam%rhs_v1_zero(1) +IF (isok) & + & CALL obj%tmp1%AXPY(x=obj%v0, scale=scale) + +scale = obj%algoParam%rhs_a1(1) * dt2 +isok = .NOT. obj%algoParam%rhs_a1_zero(1) +IF (isok) & + & CALL obj%tmp1%AXPY(x=obj%a0, scale=scale) + +scale = one +CALL obj%massMat%matVec(y=obj%rhs, x=obj%tmp1, & + addContribution=.TRUE., scale=scale) + +!---------------------------------------------------------------------------- +! make C*dt*(a*U+b*V*dt+c*A*dt*dt) +!---------------------------------------------------------------------------- + +! CALL obj%tmp1%Set(VALUE=0.0_DFP) +! +! scale = obj%algoParam%rhs_u1(2) +! isok = .NOT. obj%algoParam%rhs_u1_zero(2) +! IF (isok) & +! & CALL obj%tmp1%AXPY(x=obj%u0, scale=scale) +! +! scale = obj%algoParam%rhs_v1(2) * dt +! isok = .NOT. obj%algoParam%rhs_v1_zero(2) +! IF (isok) & +! & CALL obj%tmp1%AXPY(x=obj%v0, scale=scale) +! +! scale = obj%algoParam%rhs_a1(2) * dt2 +! isok = .NOT. obj%algoParam%rhs_a1_zero(2) +! IF (isok) & +! & CALL obj%tmp1%AXPY(x=obj%a0, scale=scale) +! +! scale = dt +! CALL obj%dampingMat%matVec(y=obj%rhs, x=obj%tmp1, & +! & addContribution=.TRUE., scale=scale) + +!---------------------------------------------------------------------------- +! make K*dt+dt*(a*U+b*V*dt+c*A*dt*dt) +!---------------------------------------------------------------------------- +CALL obj%tmp1%Set(VALUE=zero) + +scale = obj%algoParam%rhs_u1(3) +isok = .NOT. obj%algoParam%rhs_u1_zero(3) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%u0, scale=scale) + +scale = obj%algoParam%rhs_v1(3) * dt +isok = .NOT. obj%algoParam%rhs_v1_zero(3) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%v0, scale=scale) + +scale = obj%algoParam%rhs_a1(3) * dt2 +isok = .NOT. obj%algoParam%rhs_a1_zero(3) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%a0, scale=scale) + +scale = dt2 +CALL obj%stiffMat%matVec(y=obj%rhs, x=obj%tmp1, & + & addContribution=.TRUE., scale=scale) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & +& '[END] ') + +END PROCEDURE obj_AssembleRHS + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_AssembleSurfaceSource +CHARACTER(*), PARAMETER :: myName = "obj_AssembleSurfaceSource()" +REAL(DFP) :: times(1) +REAL(DFP), ALLOCATABLE :: xij(:, :), fevec(:, :), forceVec(:, :) +LOGICAL(LGT) :: isok +INTEGER(I4B) :: tnbc, ibc, idof, tmesh, imesh, id, dim, nns, iel, & + nrow, ncol, tsize, tel, con(obj%maxNNE) +INTEGER(I4B), ALLOCATABLE :: meshID(:) +CLASS(neumannBC_), POINTER :: nbc => NULL() +CLASS(AbstractMesh_), POINTER :: meshptr +TYPE(FEVariable_) :: forceVar + +isok = ALLOCATED(obj%nbc) + +IF (.NOT. isok) RETURN + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') + +CALL obj%force2%set(VALUE=zero) + +tnbc = SIZE(obj%nbc) +dim = nsd - 1_I4B +times = obj%currentTime + obj%timeStepSize(obj%currentTimeStep) + +DO ibc = 1, tnbc + nbc => obj%nbc(ibc)%ptr + isok = ASSOCIATED(nbc) + IF (.NOT. isok) CYCLE + + ! store force in tmp1 + CALL obj%tmp1%ApplyDirichletBC(dbc=nbc, times=times) + nbc => NULL() +END DO + +DO ibc = 1, tnbc + nbc => obj%nbc(ibc)%ptr + isok = ASSOCIATED(nbc) + CALL AssertError1(isok, myname, "NeumannBC is not associated") + + CALL nbc%GetParam(isSelectionByMeshID=isok) + call AssertError1(isok, myname, "currently only SelectionByMeshID is supported") + CALL nbc%GetParam(isNormal=isok) + CALL AssertError1(.NOT. isok, myname, "Normal is not supported now") + CALL nbc%GetParam(isTangent=isok) + CALL AssertError1(.NOT. isok, myname, "Tangent is not supported now") + + meshID = nbc%GetMeshID(dim=dim) + idof = nbc%GetDOFNo() + tmesh = SIZE(meshID) + + DO imesh = 1, tmesh + id = meshID(imesh) + meshptr => obj%boundaries(id)%ptr + + tel = meshptr%GetTotalElements() + + isok = meshptr%isEmpty() + IF (isok) CYCLE + + CALL obj%fedofNBC%DEALLOCATE() + CALL obj%fedofNBC%Initiate(baseContinuity=obj%baseContinuityForSpace, & + baseInterpolation=obj%baseInterpolationForSpace, & + order=obj%spaceOrder, & + mesh=meshptr, & + ipType=obj%ipTypeForSpace) + + nns = meshptr%GetMaxNNE() + CALL Reallocate(xij, 3, nns) + CALL Reallocate(fevec, nsd, nns) + CALL Reallocate(forceVec, nsd, nns) + + DO iel = 1, tel + + CALL obj%fedofNBC%GetQuadraturePoints1(quad=obj%quadForSpaceBnd, & + globalElement=iel, & + quadratureType=obj%quadTypeForSpace, & + order=2 * obj%spaceOrder, & + islocal=.TRUE.) + + CALL obj%fedofNBC%GetLocalElemShapeData(globalElement=iel, & + elemsd=obj%elemsdForSpaceBnd, & + quad=obj%quadForSpaceBnd, & + islocal=.TRUE.) + + CALL meshptr%GetNodeCoord(nodecoord=xij, nrow=nrow, ncol=ncol, & + globalElement=iel, islocal=.TRUE.) + + CALL obj%fedofNBC%GetGlobalElemShapeData(globalElement=iel, & + elemsd=obj%elemsdForSpaceBnd, & + xij=xij, islocal=.TRUE.) + + CALL obj%fedofNBC%GetConnectivity_(globalElement=iel, islocal=.TRUE., & + ans=obj%cellcon, tsize=tsize, opt="A") + con(1:tsize) = meshptr%GetGlobalNodeNumber(obj%cellcon(1:tsize)) + + CALL obj%tmp1%Get(VALUE=forceVec, globalNode=con(1:tsize), & + nrow=nrow, ncol=ncol, islocal=.TRUE., storageFMT=NODES_FMT) + + forceVar = NodalVariable(val=forceVec, rank=TypeFEVariableVector, & + vartype=TypeFEVariableSpace) + + fevec = ForceVector(test=obj%elemsdforSpaceBnd, c=forceVar, & + crank=TypeFEVariableVector) + + CALL obj%force2%Set(globalNode=con(1:tsize), VALUE=fevec, & + scale=one, addContribution=.TRUE., & + islocal=.TRUE., storageFMT=NODES_FMT) + END DO + END DO +END DO + +NULLIFY (meshptr, nbc) + +IF (ALLOCATED(meshID)) DEALLOCATE (meshID) +IF (ALLOCATED(fevec)) DEALLOCATE (fevec) +IF (ALLOCATED(forceVec)) DEALLOCATE (forceVec) +IF (ALLOCATED(xij)) DEALLOCATE (xij) + +CALL DEALLOCATE (forceVar) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & +& '[END] ') + +END PROCEDURE obj_AssembleSurfaceSource + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE AssembleMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@ConstructorMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@ConstructorMethods.F90 new file mode 100644 index 000000000..3632c9ccb --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@ConstructorMethods.F90 @@ -0,0 +1,327 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) ConstructorMethods + +USE StringUtility, ONLY: UpperCase + +USE GlobalData, ONLY: DOF_FMT, & + NONE + +USE ReallocateUtility, ONLY: Reallocate + +USE CSRMatrix_Method, ONLY: CSRMatrix_Initiate => Initiate, & + CSRMatrix_SetSparsity => SetSparsity + +USE DOF_Method, ONLY: DOF_Initiate => Initiate, & + DOF_SIZE => Size + +USE RealVector_Method, ONLY: RealVector_Initiate => Initiate, & + RealVector_Set => Set + +USE HDF5File_Class, ONLY: HDF5File_ +USE FPL +USE MatrixField_Class, ONLY: SetMatrixFieldParam +USE VectorField_Class, ONLY: SetVectorFieldParam +USE NeumannBC_Class, ONLY: NeumannBC_ + +IMPLICIT NONE + +CONTAINS + +!---------------------------------------------------------------------------- +! - obj_ImportFromToml1 +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Initiate +CHARACTER(*), PARAMETER :: myName = "obj_Initiate()" + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +CALL FPL_Init +CALL obj%param%Initiate() + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_Initiate + +!---------------------------------------------------------------------------- +! obj_ImportFromToml1 +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Deallocate +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Deallocate()" +#endif + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +obj%baseContinuityForSpace = "H1" +obj%baseInterpolationForSpace = "LAGR" +obj%verbosity = 0 + +CALL obj%param%DEALLOCATE() +CALL FPL_Finalize + +IF (ALLOCATED(obj%cellcon)) DEALLOCATE (obj%cellcon) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_Deallocate + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_InitiateDomains +CHARACTER(*), PARAMETER :: myName = "obj_InitiateDomain()" +TYPE(HDF5File_) :: meshfile +INTEGER(I4B) :: tsize, ii + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +CALL meshfile%Initiate(filename=obj%meshfilename%chars(), & + mode="READ") +CALL meshfile%OPEN() +IF (debug) CALL Display("Read meshfile") + +CALL obj%dom%Initiate(meshfile, '') +IF (debug) CALL Display("Initiate domain") + +tsize = obj%dom%GetTotalEntities(nsd - 1) +CALL Display(tsize, "tsize ::") + +ALLOCATE (obj%boundaries(tsize)) + +! NOTE: currently curve eintities is not read +! in initiate of fedomain so I load curve entities date +! independently here. +DO ii = 1, tsize + CALL Display(ii, "id ::") + ! WARN: it is necessary to allocate ptr + ALLOCATE (obj%boundaries(ii)%ptr) + CALL obj%boundaries(ii)%ptr%Initiate( & + meshfile, '', dim=nsd - 1, entities=[ii]) +END DO +IF (debug) CALL Display("Initiate Boundaries") + +CALL meshfile%DEALLOCATE() + +obj%cellmesh => obj%dom%GetMeshPointer(nsd) +obj%maxNNE = obj%cellmesh%GetMaxNNE() +obj%totalSpaceElements = obj%cellmesh%GetTotalElements() +! TODO: improve it +obj%spaceOrder = 1_I4B + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_InitiateDomains + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_InitiateFEDOF +CHARACTER(*), PARAMETER :: myName = "obj_InitiateFEDOF()" + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +CALL obj%fedof%Initiate(baseContinuity=obj%baseContinuityForSpace, & + baseInterpolation=obj%baseInterpolationForSpace, & + order=obj%spaceOrder, & + mesh=obj%cellmesh, & + ipType=obj%ipTypeForSpace) + +obj%maxCON = obj%fedof%GetMaxTotalConnectivity() +CALL Reallocate(obj%cellcon, obj%maxcon) + +! TODO: implement it +! IF (ALLOCATED(obj%nbc)) CALL InitiateFEDOF_NBC(obj) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_InitiateFEDOF + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +! SUBROUTINE InitiateFEDOF_NBC(obj) +! CLASS(UP2DFEM_), INTENT(INOUT) :: obj +! CHARACTER(*), PARAMETER :: myName = "InitiateFEDOF_NBC" +! INTEGER(I4B) :: tnbc, ibc, imesh, tmesh, id, dim +! INTEGER(I4B), ALLOCATABLE :: meshID(:) +! LOGICAL(LGT) :: isok +! CLASS(neumannBC_), POINTER :: nbc +! CLASS(AbstractMesh_), POINTER :: meshptr +! +! tnbc = SIZE(obj%nbc) +! ALLOCATE (obj%fedofNBC(tnbc)) +! +! dim = nsd - 1_I4B +! +! DO ibc = 1, tnbc +! nbc => obj%nbc(ibc)%ptr +! isok = ASSOCIATED(nbc) +! CALL AssertError1(isok, myname, "NeumannBC is not associated") +! +! CALL nbc%GetParam(isSelectionByMeshID=isok) +! CALL AssertError1(isok, myname, & +! "currently only SelectionByMeshID is supported") +! +! meshID = nbc%GetMeshID(dim=dim) +! tmesh = SIZE(meshID) +! +! DO imesh = 1, tmesh +! id = meshID(imesh) +! meshptr => obj%dom%GetMeshPointer(dim=dim, entityNum=id) +! +! isok = meshptr%isEmpty() +! IF (.NOT. isok) CYCLE +! END DO +! END DO +! +! END SUBROUTINE InitiateFEDOF_NBC + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_InitiateFields +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_InitiateFields()" +#endif + +CHARACTER(LEN=1), PARAMETER :: names(1) = ["u"] +INTEGER(I4B) :: spaceCompo(1), timeCompo(1), aint +CHARACTER(*), PARAMETER :: engineName = "NATIVE_SERIAL" + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +timeCompo(1) = 1 +spaceCompo(1) = nsd +aint = nsd * obj%maxNNE +CALL Reallocate(obj%ks, aint, aint) +CALL Reallocate(obj%ms, aint, aint) +CALL Reallocate(obj%cs, aint, aint) + +CALL SetMatrixFieldParam(param=obj%param, name="T", matrixProp="UNSYM", & + engine=engineName, spaceCompo=nsd, timeCompo=1_I4B) + +CALL obj%tanmat%Initiate(obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate tanmat") +IF (debug) CALL Display(obj%tanmat%SIZE(), 'size of tanmat:') +IF (debug) CALL Display(obj%tanmat%SHAPE(), 'shape of tanmat:') + +CALL SetMatrixFieldParam(param=obj%param, name="M", matrixProp="UNSYM", & + engine=engineName, spaceCompo=nsd, timeCompo=1_I4B) + +CALL obj%massMat%Initiate(obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate massMat") +IF (debug) CALL Display(obj%massMat%SIZE(), 'size of massMat:') +IF (debug) CALL Display(obj%massMat%SHAPE(), 'shape of massMat:') + +CALL obj%massMat%Set(VALUE=zero) + +CALL SetMatrixFieldParam(param=obj%param, name="K", matrixProp="UNSYM", & + engine=engineName, spaceCompo=nsd, timeCompo=1_I4B) + +CALL obj%stiffMat%Initiate(obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate stiffMat") +IF (debug) CALL Display(obj%stiffMat%SIZE(), 'size of stiffMat:') +IF (debug) CALL Display(obj%stiffMat%SHAPE(), 'shape of stiffMat:') + +CALL obj%stiffMat%Set(VALUE=zero) + +CALL SetVectorFieldParam(param=obj%param, name="rhs", & + engine=engineName, spaceCompo=nsd) +CALL obj%rhs%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate rhs") + +CALL SetVectorFieldParam(param=obj%param, name="sol", & + engine=engineName, spaceCompo=nsd) +CALL obj%sol%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate sol") + +CALL SetVectorFieldParam(param=obj%param, name="rhs0", & + engine=engineName, spaceCompo=nsd) +CALL obj%rhs0%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate rhs0") + +CALL SetVectorFieldParam(param=obj%param, name="u0", & + engine=engineName, spaceCompo=nsd) +CALL obj%u0%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate u0") + +CALL SetVectorFieldParam(param=obj%param, name="v0", & + engine=engineName, spaceCompo=nsd) +CALL obj%v0%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate v0") + +CALL SetVectorFieldParam(param=obj%param, name="a0", & + engine=engineName, spaceCompo=nsd) +CALL obj%a0%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate a0") + +CALL SetVectorFieldParam(param=obj%param, name="force1", & + engine=engineName, spaceCompo=nsd) +CALL obj%force1%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate force1") + +CALL SetVectorFieldParam(param=obj%param, name="force2", & + engine=engineName, spaceCompo=nsd) +CALL obj%force2%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate force2") + +CALL SetVectorFieldParam(param=obj%param, name="tmp1", & + engine=engineName, spaceCompo=nsd) +CALL obj%tmp1%Initiate(param=obj%param, fedof=obj%fedof) +IF (debug) CALL Display("Initiate tmp1") + +CALL obj%rhs%Set(VALUE=zero) +CALL obj%sol%Set(VALUE=zero) +CALL obj%rhs0%Set(VALUE=zero) +CALL obj%u0%Set(VALUE=zero) +CALL obj%v0%Set(VALUE=zero) +CALL obj%a0%Set(VALUE=zero) +CALL obj%force1%Set(VALUE=zero) +CALL obj%force2%Set(VALUE=zero) +CALL obj%tmp1%Set(VALUE=zero) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_InitiateFields + +!---------------------------------------------------------------------------- +! SetTotalDOFSpace +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE ConstructorMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@IOMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@IOMethods.F90 new file mode 100644 index 000000000..b8e61b783 --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@IOMethods.F90 @@ -0,0 +1,495 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) IOMethods + +USE TomlUtility, ONLY: GetValue, GetValue_ + +USE StringUtility, ONLY: UpperCase + +USE GlobalData, ONLY: stdout, & + CHAR_LF, & + DOF_FMT, & + NONE, & + LIS_GMRES, & + CHAR_SLASH + +USE BaseInterpolation_Method, ONLY: BaseInterpolation_ToInteger, & + BaseType_ToInteger, & + BaseType_ToChar, & + BaseInterpolation_ToChar + +USE ReallocateUtility, ONLY: Reallocate + +USE BaseType, ONLY: elem => TypeElemNameOpt, & + TypeFEVariableMatrix, TypeFEVariableConstant + +USE DirichletBC_Class, ONLY: DirichletBCImportFromToml +USE NeumannBC_Class, ONLY: NeumannBCImportFromToml +USE SolidMaterial_Class, ONLY: SolidMaterialImportFromToml +USE FEVariable_Method, ONLY: NodalVariable +USE AbstractLinSolver_Class, ONLY: AbstractLinSolverImportFromToml + +IMPLICIT NONE + +REAL(DFP), PARAMETER :: one = 1.0_DFP, zero = 0.0_DFP, minus_one = -1.0_DFP, & + half = 0.5_DFP + +CONTAINS + +!---------------------------------------------------------------------------- +! Display +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Display +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Display()" +#endif + +LOGICAL(LGT) :: isok +CHARACTER(:), ALLOCATABLE :: astr + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +CALL Display(obj%totalSpaceElements, "totalSpaceElements: ", unitno=unitno) +CALL Display(obj%totalTimeSteps, "totalTimeSteps: ", unitno=unitno) + +CALL Display(obj%baseContinuityForSpace, "baseContinuityForSpace: ", & + unitno=unitno) + +CALL Display(obj%baseInterpolationForSpace, "baseInterpolationForSpace: ", & + unitno=unitno) + +astr = BaseType_ToChar(obj%baseTypeForSpace) +CALL Display(astr, "baseTypeForSpace: ", unitno=unitno) + +astr = BaseInterpolation_TOChar(obj%ipTypeForSpace) +CALL Display(astr, "ipTypeForSpace: ", unitno=unitno) + +astr = BaseInterpolation_ToChar(obj%quadTypeForSpace) +CALL Display(astr, "quadTypeForSpace: ", unitno=unitno) + +CALL Display(obj%maxSpaceOrder, "maxSpaceOrder: ", unitno=unitno) +CALL Display(obj%timeRange, "timeRange: ", unitno=unitno) + +CALL Display(obj%spaceOrder, "spaceOrder: ", unitno=unitno) + +isok = ALLOCATED(obj%timeStepSize) +CALL Display(isok, "timeStepSize: ALLOCATED: ", unitno=unitno) +IF (isok) THEN + CALL Display(obj%timeStepSize, "timeStepSize: ", unitno=unitno) +END IF + +CALL Display(obj%result_dir%chars(), "result_dir: ", unitno=unitno) +CALL Display(obj%filename%chars(), "filename: ", unitno=unitno) + +isok = ASSOCIATED(obj%bodyForce) +CALL Display(isok, "bodyForce ASSOCIATED: ", unitno=unitno) + +isok = ASSOCIATED(obj%initialDisp) +CALL Display(isok, "initialDisp ASSOCIATED: ", unitno=unitno) + +isok = ASSOCIATED(obj%initialVel) +CALL Display(isok, "initialVel ASSOCIATED: ", unitno=unitno) + +isok = ASSOCIATED(obj%initialAcc) +CALL Display(isok, "initialAcc ASSOCIATED: ", unitno=unitno) + +CALL obj%algoParam%Display("SDAlgoParam: ", unitno=unitno) + +CALL Display("Output Controll") +CALL Display(obj%saveData, "saveData: ", unitno=unitno) +CALL Display(obj%outputFreq, "outputFreq: ", unitno=unitno) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif + +END PROCEDURE obj_Display + +!---------------------------------------------------------------------------- +! ImportFromToml +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ImportFromToml1 +CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml1()" + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START]') + +CALL ImportEssentialsFromToml(obj, table) + +CALL obj%InitiateDomains() + +CALL ImportParametersFromToml(obj, table) + +CALL ImportFunctionsFromToml(obj, table) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_ImportFromToml1 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE ImportEssentialsFromToml(obj, table) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + TYPE(toml_table), INTENT(inout) :: table + CHARACTER(*), PARAMETER :: myName = "ImportEssentialsFromToml" + INTEGER(I4B) :: origin, stat, tsize + LOGICAL(LGT) :: isok, abool + TYPE(String) :: astr + REAL(DFP), ALLOCATABLE :: temprealvec(:) + TYPE(toml_array), POINTER :: array + TYPE(toml_table), POINTER :: node + LOGICAL(LGT), ALLOCATABLE :: tempboolvec(:) + + IF (debug) CALL Display(myName//" result_dir") + + CALL GetValue(table=table, key="result_dir", VALUE=obj%result_dir, & + default_value=default_result_dir, origin=origin, stat=stat, isfound=isok) + + IF (debug) CALL Display(myName//" filename") + + CALL GetValue(table=table, key="filename", VALUE=obj%filename, & + origin=origin, stat=stat, isfound=isok, default_value=default_filename) + + CALL GetValue(table=table, key="verbosity", VALUE=obj%verbosity, & + origin=origin, stat=stat, isfound=isok, default_value=default_verbosity) + + IF (debug) CALL Display(myName//" meshfilename") + + CALL GetValue(table=table, key="meshfilename", VALUE=obj%meshfilename, & + origin=origin, stat=stat, isfound=isok, default_value="") + CALL AssertError1(isok, myname, "meshfilename not found") + + IF (debug) CALL Display(myName//" totalTimeSteps") + + CALL GetValue(table=table, key="totalTimeSteps", & + VALUE=obj%totalTimeSteps, & + default_value=0_I4B, origin=origin, stat=stat, isfound=isok) + + IF (debug) CALL Display(myName//" baseInterpolationForSpace") + + CALL GetValue(table=table, key="baseInterpolationForSpace", & + VALUE=astr, default_value=default_baseInterpolationForSpace, & + origin=origin, stat=stat, isfound=isok) + obj%baseInterpolationForSpace = UpperCase(astr%slice(1, 4)) + + CALL obj%algoParam%ImportFromToml(table=table) + + IF (debug) CALL Display(myName//" baseTypeForSpace") + + CALL GetValue(table=table, key="baseTypeForSpace", & + VALUE=astr, default_value=default_baseTypeForSpace, & + origin=origin, stat=stat, isfound=isok) + obj%baseTypeForSpace = BaseType_ToInteger(astr%chars()) + + IF (debug) CALL Display(myName//" ipTypeForSpace") + + CALL GetValue(table=table, key="ipTypeForSpace", & + VALUE=astr, default_value=default_ipTypeForSpace, & + origin=origin, stat=stat, isfound=isok) + obj%ipTypeForSpace = BaseInterpolation_ToInteger(astr%chars()) + + IF (debug) CALL Display(myName//" quadTypeForSpace") + + CALL GetValue(table=table, key="quadTypeForSpace", VALUE=astr, & + default_value=default_quadTypeForSpace, origin=origin, & + stat=stat, isfound=isok) + obj%quadTypeForSpace = BaseInterpolation_ToInteger(astr%chars()) + + IF (debug) CALL Display(myName//" timeRange") + + CALL GetValue_(table=table, key="timeRange", tsize=tsize, & + VALUE=obj%timeRange, origin=origin, stat=stat, isfound=isok) + isok = tsize .EQ. 2 + CALL AssertError1(isok, myname, "timeRange should have 2 values") + + IF (debug) CALL Display(myName//" timeStepSize") + + CALL GetValue(table=table, key="timeStepSize", VALUE=temprealvec, & + origin=origin, stat=stat, isfound=isok) + CALL AssertError1(isok, myname, "timeStepSize not found") + + CALL Reallocate(obj%timeStepSize, obj%totalTimeSteps) + + abool = SIZE(temprealvec) .EQ. 1 + IF (abool) THEN + obj%timeStepSize = temprealvec(1) + ELSE + isok = SIZE(temprealvec) .EQ. obj%totalTimeSteps + CALL AssertError1(isok, myname, "timeStepSize should have "// & + "totalTimeElements values") + obj%timeStepSize(:) = temprealvec(1:obj%totalTimeSteps) + END IF + + IF (debug) CALL Display(myName//" saveData") + + obj%saveData = .TRUE. + CALL toml_get(table, "saveData", array, origin=origin, stat=stat) + CALL toml_get(array, tempboolvec, origin=origin, stat=stat) + + abool = SIZE(tempboolvec) .EQ. 1 + IF (abool) THEN + obj%saveData = tempboolvec(1) + ELSE + isok = SIZE(tempboolvec) .EQ. 4 + CALL AssertError1(isok, myname, "saveData should have 4 values") + obj%saveData = tempboolvec + END IF + array => NULL() + DEALLOCATE (tempboolvec) + + IF (debug) CALL Display(myName//" outputFreq") + + CALL GetValue(table=table, key="outputFreq", VALUE=obj%outputFreq, & + default_value=1_I4B, origin=origin, stat=stat) + + IF (debug) CALL Display(myName//" Linsolver") + CALL toml_get(table, "linsolver", node, origin=origin, & + requested=.FALSE., stat=stat) + isok = ASSOCIATED(node) + CALL asserterror1(isok, myname, "linsolver is not found") + CALL AbstractLinSolverImportFromToml(obj%linsolver, node) + + node => NULL() + +END SUBROUTINE ImportEssentialsFromToml + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE ImportParametersFromToml(obj, table) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + TYPE(toml_table), INTENT(inout) :: table + CHARACTER(*), PARAMETER :: myName = "ImportParametersFromToml" + INTEGER(I4B) :: origin, stat, tsize + LOGICAL(LGT) :: isok + TYPE(toml_table), POINTER :: node + TYPE(toml_array), POINTER :: array + TYPE(UserFunction_), POINTER :: func => NULL() + REAL(DFP) :: mat(3, 3) + + ! TODO: need to improve this method + ! for example I want to set the differenct values for + ! specific domains + IF (debug) CALL Display(myName//" Solid material") + node => NULL() + CALL toml_get(table, "solidMaterial", array, origin=origin, & + requested=.FALSE., stat=stat) + isok = ASSOCIATED(array) + CALL asserterror1(isok, myname, "solidMaterial not found") + tsize = toml_len(array) + isok = tsize .EQ. 1 + CALL asserterror1(isok, myname, & + "currently only homogeneous material is supported") + ALLOCATE (obj%solidMaterial(tsize)) + IF (debug) CALL Display(tsize, "number of solidMaterial :: ") + CALL SolidMaterialImportFromToml(obj=obj%solidMaterial, & + tomlName="solidMaterial", table=table) + array => NULL() + + func => obj%solidMaterial(1)%ptr%GetMaterialPointer("massDensity") + isok = ASSOCIATED(func) + CALL AssertError1(isok, myname, "massDensity not found") + CALL func%GetFEVariable(obj%density) + func => NULL() + + func => obj%solidMaterial(1)%ptr%GetMaterialPointer("rayleigh_alpha") + isok = ASSOCIATED(func) + CALL AssertError1(isok, myname, "rayleigh_alpha not found") + CALL func%GetFEVariable(obj%rayleighAlpha) + func => NULL() + + func => obj%solidMaterial(1)%ptr%GetMaterialPointer("rayleigh_beta") + isok = ASSOCIATED(func) + CALL AssertError1(isok, myname, "rayleigh_beta not found") + CALL func%GetFEVariable(obj%rayleighBeta) + func => NULL() + + CALL obj%solidMaterial(1)%ptr%stressStrainModel%GetC(mat) + obj%cijkl = NodalVariable(mat, & + & TypeFEVariableMatrix, & + & TypeFEVariableConstant) + +END SUBROUTINE ImportParametersFromToml + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +SUBROUTINE ImportFunctionsFromToml(obj, table) + CLASS(UP2DFEM_), INTENT(INOUT) :: obj + TYPE(toml_table), INTENT(inout) :: table + CHARACTER(*), PARAMETER :: myName = "ImportFunctionsFromToml" + INTEGER(I4B) :: origin, stat, tsize + LOGICAL(LGT) :: isok + TYPE(String) :: astr + TYPE(toml_table), POINTER :: node + TYPE(toml_array), POINTER :: array + + astr = "bodyForce" + IF (debug) CALL Display(myName//astr%chars()) + node => NULL() + CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & + stat=stat) + isok = ASSOCIATED(node) + IF (isok) THEN + ALLOCATE (obj%bodyForce) + CALL obj%bodyForce%ImportFromToml(table=node) + END IF + node => NULL() + + astr = "initialAcc" + IF (debug) CALL Display(myName//astr%chars()) + node => NULL() + CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & + stat=stat) + isok = ASSOCIATED(node) + IF (isok) THEN + ALLOCATE (obj%initialAcc) + CALL obj%initialAcc%ImportFromToml(table=node) + END IF + node => NULL() + + astr = "initialVel" + IF (debug) CALL Display(myName//astr%chars()) + node => NULL() + CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & + stat=stat) + isok = ASSOCIATED(node) + IF (isok) THEN + ALLOCATE (obj%initialVel) + CALL obj%initialVel%ImportFromToml(table=node) + END IF + node => NULL() + + astr = "initialDisp" + IF (debug) CALL Display(myName//astr%chars()) + node => NULL() + CALL toml_get(table, astr%chars(), node, origin=origin, requested=.FALSE., & + stat=stat) + isok = ASSOCIATED(node) + IF (isok) THEN + ALLOCATE (obj%initialDisp) + CALL obj%initialDisp%ImportFromToml(table=node) + END IF + node => NULL() + + astr = "dirichletBC" + IF (debug) CALL Display(myName//astr%chars()) + node => NULL() + CALL toml_get(table, astr%chars(), array, origin=origin, requested=.FALSE., & + stat=stat) + isok = ASSOCIATED(array) + IF (isok) THEN + tsize = toml_len(array) + ALLOCATE (obj%dbc(tsize)) + IF (debug) CALL Display(tsize, "number of DirichletBC :: ") + CALL DirichletBCImportFromToml(obj=obj%dbc, tomlName="dirichletBC", & + table=table, dom=obj%dom) + array => NULL() + ELSE + IF (debug) CALL Display("No dirichletBC") + END IF + + astr = "neumannBC" + IF (debug) CALL Display(myName//astr%chars()) + node => NULL() + CALL toml_get(table, astr%chars(), array, origin=origin, requested=.FALSE., & + stat=stat) + isok = ASSOCIATED(array) + IF (isok) THEN + tsize = toml_len(array) + ALLOCATE (obj%nbc(tsize)) + IF (debug) CALL Display(tsize, "number of NeumannBC :: ") + CALL NeumannBCImportFromToml(obj=obj%nbc, tomlName="neumannBC", & + table=table, dom=obj%dom) + array => NULL() + ELSE + IF (debug) CALL Display("No neumannBC") + END IF + +END SUBROUTINE ImportFunctionsFromToml + +!---------------------------------------------------------------------------- +! ImportFromToml +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ImportFromToml2 +! internal variables +CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml2()" +TYPE(toml_table), ALLOCATABLE :: table +TYPE(toml_table), POINTER :: node +INTEGER(I4B) :: origin, stat +LOGICAL(LGT) :: isok + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START]') +#endif + +CALL GetValue(table=table, afile=afile, filename=filename) + +isok = ALLOCATED(table) +CALL AssertError1(isok, myname, "table is not allocated from GetValue") + +node => NULL() +CALL toml_get(table, tomlName, node, origin=origin, requested=.FALSE., & + stat=stat) + +IF (.NOT. ASSOCIATED(node)) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + '[INTERNAL ERROR] :: following error occured while reading '// & + 'the toml file :: cannot find '//tomlName//" table in config.") +END IF + +CALL obj%ImportFromToml(table=node) + +#ifdef DEBUG_VER +IF (PRESENT(printToml)) THEN + CALL Display(toml_serialize(node), myname//" Domain toml config: "// & + CHAR_LF, unitno=stdout) +END IF +#endif + +node => NULL() + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END]') +#endif + +END PROCEDURE obj_ImportFromToml2 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE IOMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@RunMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@RunMethods.F90 new file mode 100644 index 000000000..e72810f69 --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@RunMethods.F90 @@ -0,0 +1,57 @@ +SUBMODULE(UP2DFEM_Class) RunMethods + +IMPLICIT NONE + +CONTAINS + +!---------------------------------------------------------------------------- +! Solve +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Run +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Run()" +#endif +INTEGER(I4B) :: itime + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') + +IF (debug) CALL Display("Simulation Starts") + +CALL obj%SetInitialAcceleration() +CALL obj%SetInitialVelocity() +CALL obj%SetInitialDisplacement() + +CALL obj%AssembleTanmat() +CALL obj%WriteData() + +DO itime = 1, obj%totalTimeSteps + + CALL Display(obj%currentTime, myname//" current time: ") + CALL obj%AssembleTanmat() + ! CALL obj%AssembleBodySource() + CALL obj%AssembleSurfaceSource() + ! CALL obj%AssemblePointSource() + CALL obj%AssembleRHS() + CALL obj%ApplyDirichletBC() + CALL obj%Solve() + CALL obj%Update() + CALL obj%WriteData() + +END DO + +IF (debug) CALL Display("Simulation is Finished") + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & +& '[END]') + +END PROCEDURE obj_Run + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE RunMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@SetMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@SetMethods.F90 new file mode 100644 index 000000000..1a7c3df81 --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@SetMethods.F90 @@ -0,0 +1,154 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) SetMethods + +IMPLICIT NONE + +REAL(DFP), PARAMETER :: one = 1.0_DFP, zero = 0.0_DFP, minus_one = -1.0_DFP, & + half = 0.5_DFP + +CONTAINS + +!---------------------------------------------------------------------------- +! obj_Set +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Set +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Set()" +#endif +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +obj%currentTimeStep = 1 +obj%currentTime = obj%timeRange(1) + +CALL obj%InitiateFEDOF() +CALL obj%InitiateFields() + +isok = ALLOCATED(obj%dbc) + +IF (isok) THEN + obj%dbcIndx = obj%sol%GetNodeLoc(dbc=obj%dbc) +END IF + +CALL obj%linsolver%set(obj%tanmat) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_Set + +!---------------------------------------------------------------------------- +! setInitialAcceleration +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_SetInitialAcceleration +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_SetInitialAcceleration()" +#endif + +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +isok = ASSOCIATED(obj%initialAcc) +IF (.NOT. isok) THEN + CALL obj%a0%set(VALUE=0.0_DFP) + RETURN +END IF + +! WARN: currently initial acceleration is directly set with +! nodal value. It does not support all interpolation +CALL obj%a0%SetByFunction(obj%initialAcc) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_SetInitialAcceleration + +!---------------------------------------------------------------------------- +! SetInitialVelocity +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_SetInitialVelocity +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_SetInitialVelocity()" +#endif + +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +isok = ASSOCIATED(obj%initialVel) +IF (.NOT. isok) THEN + CALL obj%v0%set(VALUE=0.0_DFP) + RETURN +END IF + +! WARN: currently initial value is directly set with +! nodal value. It does not support all interpolation +CALL obj%v0%SetByFunction(obj%initialVel) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_SetInitialVelocity + +!---------------------------------------------------------------------------- +! SetInitialVelocity +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_SetInitialDisplacement +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_SetInitialDisplacement()" +#endif + +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +isok = ASSOCIATED(obj%initialDisp) +IF (.NOT. isok) THEN + CALL obj%u0%set(VALUE=0.0_DFP) + RETURN +END IF + +! WARN: currently initial value is directly set with +! nodal value. It does not support all interpolation +CALL obj%u0%SetByFunction(obj%initialDisp) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_SetInitialDisplacement + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE SetMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@SolveMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@SolveMethods.F90 new file mode 100644 index 000000000..aea552f14 --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@SolveMethods.F90 @@ -0,0 +1,51 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) SolveMethods + +IMPLICIT NONE + +CONTAINS + +!---------------------------------------------------------------------------- +! Solve +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Solve +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Solve()" +#endif + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +CALL obj%linsolver%Solve(sol=obj%sol, rhs=obj%rhs) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_Solve + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE SolveMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@UpdateMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@UpdateMethods.F90 new file mode 100644 index 000000000..84cebcbcf --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@UpdateMethods.F90 @@ -0,0 +1,152 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) UpdateMethods + +IMPLICIT NONE + +REAL(DFP), PARAMETER :: one = 1.0_DFP, zero = 0.0_DFP, minus_one = -1.0_DFP, & + half = 0.5_DFP + +CONTAINS + +!---------------------------------------------------------------------------- +! Update +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Update +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Update()" +#endif + +REAL(DFP) :: scale, dt, dts +LOGICAL(LGT) :: isok + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') + +dt = obj%timeStepSize(obj%currentTimeStep) +dts = dt * dt +obj%currentTime = obj%currentTime + dt +obj%currentTimeStep = obj%currentTimeStep + 1 + +CALL obj%force1%COPY(obj%force2) + +!---------------------------------------------------------------------------- +! Update displacement +!---------------------------------------------------------------------------- + +CALL obj%rhs%Set(VALUE=zero) + +scale = obj%algoParam%dis(1) +isok = .NOT. obj%algoParam%dis_zero(1) +IF (isok) & + CALL obj%rhs%AXPY(x=obj%u0, scale=scale) + +scale = obj%algoParam%dis(2) * dt +isok = .NOT. obj%algoParam%dis_zero(2) +IF (isok) & + CALL obj%rhs%AXPY(x=obj%v0, scale=scale) + +scale = obj%algoParam%dis(3) * dts +isok = .NOT. obj%algoParam%dis_zero(3) +IF (isok) & + CALL obj%rhs%AXPY(x=obj%a0, scale=scale) + +scale = obj%algoParam%dis(4) +isok = .NOT. obj%algoParam%dis_zero(4) +IF (isok) & + CALL obj%rhs%AXPY(x=obj%sol, scale=scale) + +!---------------------------------------------------------------------------- +! Update v0 +!---------------------------------------------------------------------------- + +CALL obj%tmp1%Set(VALUE=zero) + +scale = obj%algoParam%vel(1) / dt +isok = .NOT. obj%algoParam%vel_zero(1) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%u0, scale=scale) + +scale = obj%algoParam%vel(2) +isok = .NOT. obj%algoParam%vel_zero(2) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%v0, scale=scale) + +scale = obj%algoParam%vel(3) * dt +isok = .NOT. obj%algoParam%vel_zero(3) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%a0, scale=scale) + +scale = obj%algoParam%vel(4) / dt +isok = .NOT. obj%algoParam%vel_zero(4) +IF (isok) & + CALL obj%tmp1%AXPY(x=obj%sol, scale=scale) + +!---------------------------------------------------------------------------- +! Update a0 +!---------------------------------------------------------------------------- + +CALL obj%force2%Set(VALUE=zero) + +scale = obj%algoParam%acc(1) / dts +isok = .NOT. obj%algoParam%acc_zero(1) +IF (isok) & + CALL obj%force2%AXPY(x=obj%u0, scale=scale) + +scale = obj%algoParam%acc(2) / dt +isok = .NOT. obj%algoParam%acc_zero(2) +IF (isok) & + CALL obj%force2%AXPY(x=obj%v0, scale=scale) + +scale = obj%algoParam%acc(3) +isok = .NOT. obj%algoParam%acc_zero(3) +IF (isok) & + CALL obj%force2%AXPY(x=obj%a0, scale=scale) + +scale = obj%algoParam%acc(4) / dts +isok = .NOT. obj%algoParam%acc_zero(4) +IF (isok) & + CALL obj%force2%AXPY(x=obj%sol, scale=scale) + +!---------------------------------------------------------------------------- +! Update u0, v0, a0, and force1 +!---------------------------------------------------------------------------- + +CALL obj%u0%Copy(obj%rhs) +CALL obj%v0%Copy(obj%tmp1) +CALL obj%a0%Copy(obj%force2) +CALL obj%force2%Set(VALUE=zero) +CALL obj%rhs%Set(VALUE=zero) + +CALL obj%sol%Set(VALUE=zero) + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') + +END PROCEDURE obj_Update + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE UpdateMethods diff --git a/src/submodules/UP2DFEM/src/UP2DFEM_Class@WriteDataMethods.F90 b/src/submodules/UP2DFEM/src/UP2DFEM_Class@WriteDataMethods.F90 new file mode 100644 index 000000000..198a68c13 --- /dev/null +++ b/src/submodules/UP2DFEM/src/UP2DFEM_Class@WriteDataMethods.F90 @@ -0,0 +1,91 @@ +! This program is a part of EASIFEM library +! Expandable And Scalable Infrastructure for Finite Element Methods +! htttps://www.easifem.com +! Vikas Sharma, Ph.D., vickysharma0812@gmail.com +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see +! + +SUBMODULE(UP2DFEM_Class) WriteDataMethods + +USE VTKFile_Class, ONLY: VTKFile_, VTK_BINARY_APPENDED, & + VTK_UnstructuredGrid + +USE AbstractNodeField_Class, ONLY: AbstractNodeFieldPointer_ + +IMPLICIT NONE + +REAL(DFP), PARAMETER :: one = 1.0_DFP, zero = 0.0_DFP, minus_one = -1.0_DFP, & + half = 0.5_DFP + +CONTAINS + +!---------------------------------------------------------------------------- +! WriteData +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_WriteData +CHARACTER(*), PARAMETER :: myName = "obj_WriteData()" +TYPE(VTKFile_) :: avtk +TYPE(String) :: path, filename, suffix +CHARACTER(*), PARAMETER :: ext = ".vtu" +LOGICAL(LGT) :: isOK + +IF (debug) CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START]') + +isOk = MOD(obj%currentTimeStep - 1, obj%OutputFreq) .NE. 0_I4B +IF (isOk) RETURN + +suffix = tostring(obj%currentTimeStep - 1) +path = obj%result_dir + +isok = obj%saveData(1) +IF (isok) THEN + filename = path%chars()//obj%filename//"_disp_"//suffix//ext + CALL avtk%InitiateVTKFile(filename=filename%chars(), mode="NEW", & + DataFormat=VTK_BINARY_APPENDED, DataStructureType=VTK_UnStructuredGrid) + CALL obj%u0%WriteData(vtk=avtk) + CALL avtk%DEALLOCATE() +END IF + +isok = obj%saveData(2) +IF (isok) THEN + filename = path%chars()//obj%filename//"_velo_"//suffix//ext + CALL avtk%InitiateVTKFile(filename=filename%chars(), mode="NEW", & + DataFormat=VTK_BINARY_APPENDED, DataStructureType=VTK_UnStructuredGrid) + CALL obj%v0%WriteData(vtk=avtk) + CALL avtk%DEALLOCATE() +END IF + +isok = obj%saveData(3) +IF (isok) THEN + filename = path%chars()//obj%filename//"_acc_"//suffix//ext + CALL avtk%InitiateVTKFile(filename=filename%chars(), mode="NEW", & + DataFormat=VTK_BINARY_APPENDED, DataStructureType=VTK_UnStructuredGrid) + CALL obj%a0%WriteData(vtk=avtk) + CALL avtk%DEALLOCATE() +END IF + +IF (debug) CALL e%RaiseInformation(modName//'::'// & + myName//' - '//'[END]') +END PROCEDURE obj_WriteData + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE WriteDataMethods From cb9db069e8a05c24349dfc6711cc2d064ab1ff8d Mon Sep 17 00:00:00 2001 From: shion Date: Tue, 4 Nov 2025 12:26:42 +0900 Subject: [PATCH 2/2] backup fix --- .../LinSolver/src/LIS_SOLVE_PRECOND.F90 | 67 ------------ .../src/LinSolver_Class@SolveMethods.F90 | 100 ------------------ 2 files changed, 167 deletions(-) delete mode 100644 src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 diff --git a/src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 b/src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 deleted file mode 100644 index c80b5e7d6..000000000 --- a/src/submodules/LinSolver/src/LIS_SOLVE_PRECOND.F90 +++ /dev/null @@ -1,67 +0,0 @@ -SUBROUTINE _SUBROUTINE_NAME(obj, sol, rhs, precond) - CLASS(LinSolver_), TARGET, INTENT(INOUT) :: obj - TYPE(CSRMatrix_), INTENT(IN) :: precond - REAL(DFP), INTENT(INOUT) :: sol(:) - REAL(DFP), INTENT(INOUT) :: rhs(:) - - ! Internal variables - CHARACTER(*), PARAMETER :: myName = _MY_NAME - INTEGER(I4B) :: n - REAL(DFP), ALLOCATABLE :: diag(:) - CLASS(AbstractMatrixField_), POINTER :: amat - LOGICAL(LGT) :: isok - -#ifdef DEBUG_VER - CALL e%RaiseInformation(modName//'::'//myName//' - '// & - '[START] ') -#endif - - obj%IPAR(1) = 0 - obj%FPAR(11) = 0.0_DFP - CALL obj%GetParam(globalNumRow=n, amat=amat) - obj%IPAR(7) = 1 - - DO - - CALL _LIS_NAME(n, rhs, sol, obj%IPAR, obj%FPAR, obj%W) - - IF (obj%IPAR(1) .GT. 0) THEN - - CALL PERFORM_TASK_PRECOND(amat, y=obj%W(obj%IPAR(9):obj%IPAR(9) + n - 1), & - x=obj%W(obj%IPAR(8):obj%IPAR(8) + n - 1), & - precond=precond, & - ierr=obj%IPAR(1)) - - ELSE IF (obj%IPAR(1) .LT. 0) THEN - - CALL CHECKERROR(IPAR=obj%IPAR, FPAR=obj%FPAR, myName=myName) - EXIT - - ELSE IF (obj%IPAR(1) .EQ. 0) THEN - - CALL obj%SetParam(ierr=obj%ipar(1), iter=obj%ipar(7)) - CALL DisplayConvergence(myName, obj%ipar(7), obj%FPAR) - EXIT - - END IF - - END DO - - ! Initial residual/error norm - - CALL obj%SetParam(error0=obj%fpar(3), tol=obj%fpar(4), & - error=obj%fpar(6), normRes=obj%fpar(5)) - -END SUBROUTINE _SUBROUTINE_NAME - -#ifdef _SUBROUTINE_NAME -#undef _SUBROUTINE_NAME -#endif - -#ifdef _LIS_NAME -#undef _LIS_NAME -#endif - -#ifdef _MY_NAME -#undef _MY_NAME -#endif diff --git a/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 b/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 index ea4cb3d43..2901aefe3 100644 --- a/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 +++ b/src/submodules/LinSolver/src/LinSolver_Class@SolveMethods.F90 @@ -82,59 +82,6 @@ SUBROUTINE PERFORM_TASK(amat, y, x, ierr) END SUBROUTINE PERFORM_TASK -!---------------------------------------------------------------------------- -! -!---------------------------------------------------------------------------- - -SUBROUTINE PERFORM_TASK_PRECOND(amat, y, x, precond, ierr) - - ! intent of dummy variables - CLASS(AbstractMatrixField_), INTENT(INOUT) :: amat - TYPE(CSRMatrix_), INTENT(IN) :: precond - REAL(DFP), INTENT(INOUT) :: y(:) - REAL(DFP), INTENT(IN) :: x(:) - INTEGER(I4B), INTENT(IN) :: ierr - -#ifdef DEBUG_VER - CHARACTER(*), PARAMETER :: myName = "PERFORM_TASK()" -#endif - -#ifdef DEBUG_VER - CALL e%RaiseInformation(modName//'::'//myName//' - '// & - '[START] ') -#endif - - SELECT CASE (ierr) - - CASE (1) - - CALL amat%Matvec(y=y, x=x, isTranspose=.FALSE.) - - CASE (2) - - CALL amat%Matvec(y=y, x=x, isTranspose=.TRUE.) - - CASE (3, 5) - - ! LEFT/RIGHT PRECONDITIONER SOLVER - ! The preconditioners are given precond - CALL Matvec(obj=precond, x=x, y=y, isTranspose=.FALSE.) - - CASE (4, 6) - - ! LEFT/RIGHT PRECONDITIONER SOLVER - ! The preconditioners are given precond - CALL Matvec(obj=precond, x=x, y=y, isTranspose=.TRUE.) - - END SELECT - -#ifdef DEBUG_VER - CALL e%RaiseInformation(modName//'::'//myName//' - '// & - '[END] ') -#endif - -END SUBROUTINE PERFORM_TASK_PRECOND - !---------------------------------------------------------------------------- ! CHECKERR !---------------------------------------------------------------------------- @@ -385,43 +332,6 @@ END SUBROUTINE DisplayConvergence END PROCEDURE obj_Solve -!---------------------------------------------------------------------------- -! -!---------------------------------------------------------------------------- - -MODULE PROCEDURE obj_Solve_precond -CHARACTER(*), PARAMETER :: myName = "obj_Solve()" -REAL(DFP), POINTER :: rhsvar(:), solvar(:) -INTEGER(I4B) :: info -INTEGER(I4B) :: solverName -LOGICAL(LGT) :: isok -CLASS(AbstractMatrixField_), POINTER :: amat - -CALL obj%GetParam(isInitiated=isok, solverName=solverName, amat=amat) - -CALL AssertError1(isok, myname, 'Linear solver is not initiated!') - -isok = ASSOCIATED(amat) -CALL AssertError1(isok, myname, 'Amat is not associated') - -SELECT CASE (solverName) - -CASE (TypeSolverNameOpt%GMRES) - - rhsvar => rhs%GetPointer() - solvar => sol%GetPointer() - CALL LS_SOLVE_GMRES_precond(obj=obj, sol=solvar, rhs=rhsvar, & - precond=precond) - NULLIFY (rhsvar, solvar) - -CASE DEFAULT - - CALL AssertError1(.FALSE., myName, 'No case found for linear solver') - RETURN -END SELECT - -END PROCEDURE obj_Solve_precond - !---------------------------------------------------------------------------- ! LS_SOLVE_CG !---------------------------------------------------------------------------- @@ -495,16 +405,6 @@ END SUBROUTINE DisplayConvergence #include "./LIS_SOLVE.F90" -!---------------------------------------------------------------------------- -! -!---------------------------------------------------------------------------- - -#define _SUBROUTINE_NAME LS_SOLVE_GMRES_PRECOND -#define _LIS_NAME GMRES -#define _MY_NAME "LS_SOLVE_GMRES_PRECOND" - -#include "./LIS_SOLVE_PRECOND.F90" - !---------------------------------------------------------------------------- ! LS_SOLVE_FGMRES !----------------------------------------------------------------------------