From 9cb23c28ca09f311c7f239affa9ad2afa6952491 Mon Sep 17 00:00:00 2001 From: shion <106575883+shishiousan@users.noreply.github.com> Date: Fri, 5 Jan 2024 14:19:28 +0900 Subject: [PATCH 1/5] Updates in AcousticKernel - working in progress --- src/modules/AbstractKernel/CMakeLists.txt | 1 + .../src/AbstractKernel_Class.F90 | 30 ++++ .../KernelAssembleDiffusionMatrix_Method.F90 | 162 ++++++++++++++++++ ...ractKernel_Class@AssembleTanmatMethods.F90 | 35 ++++ 4 files changed, 228 insertions(+) create mode 100644 src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 diff --git a/src/modules/AbstractKernel/CMakeLists.txt b/src/modules/AbstractKernel/CMakeLists.txt index 06156f03c..5b7b157f0 100644 --- a/src/modules/AbstractKernel/CMakeLists.txt +++ b/src/modules/AbstractKernel/CMakeLists.txt @@ -22,6 +22,7 @@ target_sources( ${src_path}/AbstractKernel_Class.F90 ${src_path}/KernelAssembleMassMatrix_Method.F90 ${src_path}/KernelAssembleStiffnessMatrix_Method.F90 + ${src_path}/KernelAssembleDiffusionMatrix_Method.F90 ${src_path}/KernelAssembleDampingMatrix_Method.F90 ${src_path}/KernelAssembleElastoDynaMatrix_Method.F90 ${src_path}/KernelAssembleSTElastoDynaMatrix_Method.F90 diff --git a/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 b/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 index 45cc784d2..621d443ea 100644 --- a/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 +++ b/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 @@ -415,6 +415,8 @@ MODULE AbstractKernel_Class !! List of space-time scalar fields CLASS(MatrixField_), POINTER :: stiffnessMat => NULL() !! Global Stiffness matrix + CLASS(MatrixField_), POINTER :: diffusionMat => NULL() + !! Global Stiffness matrix CLASS(MatrixField_), POINTER :: massMat => NULL() !! Global mass matrix CLASS(MatrixField_), POINTER :: dampingMat => NULL() @@ -427,6 +429,12 @@ MODULE AbstractKernel_Class !! Vector field for nodal acceleration CLASS(VectorField_), POINTER :: nodeCoord => NULL() !! Vector field for nodal coordinates + CLASS(ScalarField_), POINTER :: pressure => NULL() + !! scalar field for nodal pressure + CLASS(ScalarField_), POINTER :: p_velocity => NULL() + !! scalar field for nodal pressure + CLASS(ScalarField_), POINTER :: p_acceleration => NULL() + !! scalar field for nodal pressure TYPE(VectorMeshFieldPointer_), ALLOCATABLE :: solidMechData(:) !! Constitutive data for solid materials TYPE(AbstractScalarMeshFieldPointer_), ALLOCATABLE :: massDensity(:) @@ -454,6 +462,12 @@ MODULE AbstractKernel_Class TYPE(AbstractVectorMeshFieldPointer_), ALLOCATABLE :: strain(:) !! Strain tensor !! This will be a tensor mesh field + ! TYPE(AbstractScalarMeshFieldPointer_), ALLOCATABLE :: phase_velocity(:) + !! phase_velocity + !! This will be a scalar mesh field + TYPE(AbstractScalarMeshFieldPointer_), ALLOCATABLE :: scalarCoefficient(:) + !! it can be phase velocity or coefficient of permiabillity for isotropic medium + !! this will be a scalar mesh field CLASS(UserFunction_), POINTER :: bodySourceFunc => NULL() !! body force function @@ -678,6 +692,8 @@ MODULE AbstractKernel_Class PROCEDURE, PUBLIC, PASS(obj) :: AssembleMassMat => obj_AssembleMassMat PROCEDURE, PUBLIC, PASS(obj) :: AssembleStiffnessMat => & & obj_AssembleStiffnessMat + PROCEDURE, PUBLIC, PASS(obj) :: AssembleDiffusionMat => & + & obj_AssembleDiffusionMat PROCEDURE, PUBLIC, PASS(obj) :: AssembleDampingMat => & & obj_AssembleDampingMat PROCEDURE, PUBLIC, PASS(obj) :: AssembleNitscheMat => & @@ -1870,6 +1886,20 @@ MODULE SUBROUTINE obj_AssembleStiffnessMat(obj) END SUBROUTINE obj_AssembleStiffnessMat END INTERFACE +!---------------------------------------------------------------------------- +! AssembleDiffusionMat@AssembleTanmatMethods +!---------------------------------------------------------------------------- + +!> authors: Shion Shimizu +! date: 2024-01-04 +! summary: This subroutine assembles the diffusion matrix of the system + +INTERFACE + MODULE SUBROUTINE obj_AssembleDiffusionMat(obj) + CLASS(AbstractKernel_), INTENT(INOUT) :: obj + END SUBROUTINE obj_AssembleDiffusionMat +END INTERFACE + !---------------------------------------------------------------------------- ! AssembleDampingMat@AssembleTanmatMethods !---------------------------------------------------------------------------- diff --git a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 new file mode 100644 index 000000000..c9fb5ee6f --- /dev/null +++ b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 @@ -0,0 +1,162 @@ +! 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 +! + +MODULE KernelAssembleDiffusionMatrix_Method +USE GlobalData +USE Field +USE BaseMethod +USE BaseType +USE FieldFactory +USE Domain_Class +USE Mesh_Class +USE FiniteElement_Class +USE ExceptionHandler_Class, ONLY: e +USE AbstractKernelParam, ONLY: KernelProblemType +IMPLICIT NONE +PRIVATE + +PUBLIC :: KernelAssembleDiffusionMatrix + +CHARACTER(*), PARAMETER :: modName = "KernelAssembleDiffusionMatrix_Method" + +INTERFACE KernelAssembleDiffusionMatrix + MODULE PROCEDURE KernelAssembleIsoDiffMat +END INTERFACE KernelAssembleDiffusionMatrix + +! TODO: add some other assemble method e.g. anisotropic + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +CONTAINS + +!---------------------------------------------------------------------------- +! KernelAssembleIsotropicDiffusionMatrix +!---------------------------------------------------------------------------- + +SUBROUTINE KernelAssembleIsoDiffMat(mat, coefficient, dom, cellFE, & + & linCellFE, spaceElemSD, linSpaceElemSD, reset) + CLASS(MatrixField_), INTENT(INOUT) :: mat + CLASS(AbstractScalarMeshFieldPointer_), INTENT(INOUT) :: coefficient(:) + CLASS(Domain_), INTENT(INOUT) :: dom + TYPE(FiniteElementPointer_), INTENT(INOUT) :: cellFE(:) + TYPE(FiniteElementPointer_), INTENT(INOUT) :: linCellFE(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: spaceElemSD(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: linSpaceElemSD(:) + LOGICAL(LGT), INTENT(IN) :: reset + + ! internal variables + CHARACTER(*), PARAMETER :: myName = "KernelAssembleIsoDiffMat()" + INTEGER(I4B) :: id, tmesh, nsd, nns, tdof, iel + TYPE(ElemShapeData_) :: elemsd, linElemSD + TYPE(FEVariable_) :: coefVar + CLASS(Mesh_), POINTER :: meshptr + CLASS(FiniteElement_), POINTER :: spaceFE, linSpaceFE + CLASS(ReferenceElement_), POINTER :: refelem + CLASS(AbstractScalarMeshField_), POINTER :: coefficientField + REAL(DFP), ALLOCATABLE :: Mmat(:, :), xij(:, :) + INTEGER(I4B), ALLOCATABLE :: nptrs(:) + LOGICAL(LGT) :: isok, problem + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + + IF (reset) CALL mat%Set(VALUE=0.0_DFP) + nsd = dom%GetNSD() + tmesh = dom%GetTotalMesh(dim=nsd) + tdof = 1_I4B + + NULLIFY (meshptr, spaceFE, linSpaceFE, coefficientField, refelem) + + DO id = 1, tmesh + meshptr => dom%GetMeshPointer(dim=nsd, entityNum=id) + + isok = ASSOCIATED(meshptr) + IF (.NOT. isok) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: Null mesh found at id = '//tostring(id)) + RETURN + END IF + + problem = meshptr%isEmpty() + IF (problem) CYCLE + + spaceFE => cellFE(id)%ptr + linSpaceFE => linCellFE(id)%ptr + + elemsd = spaceElemSD(id) + linElemSD = linSpaceElemSD(id) + + refelem => meshptr%GetRefElemPointer() + nns = (.NNE.refelem) + + CALL Reallocate(xij, nsd, nns) + CALL Reallocate(Mmat, tdof * nns, tdof * nns) + + coefficientField => coefficient(id)%ptr + + isok = ASSOCIATED(coefficientField) + IF (.NOT. isok) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: coefficient('//tostring(id)//') is NULL.') + RETURN + END IF + + DO iel = meshptr%minElemNum, meshptr%maxElemNum + isok = meshptr%isElementPresent(iel) + IF (.NOT. isok) CYCLE + + CALL coefficientField%Get(fevar=coefVar, globalElement=iel) + + nptrs = meshptr%GetConnectivity(iel) + + CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) + + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + & geoElemSD=linElemSD) + + Mmat = DiffusionMatrix(test=elemsd, trial=elemsd, k=coefVar, krank=coefVar) + + CALL mat%Set(globalNode=nptrs, VALUE=Mmat, & + & storageFMT=DOF_FMT, scale=1.0_DFP, addContribution=.TRUE.) + + END DO + END DO + + NULLIFY (meshptr, spaceFE, coefficientField, refelem, linSpaceFE) + + IF (ALLOCATED(Mmat)) DEALLOCATE (Mmat) + IF (ALLOCATED(xij)) DEALLOCATE (xij) + IF (ALLOCATED(nptrs)) DEALLOCATE (nptrs) + CALL DEALLOCATE (elemsd) + CALL DEALLOCATE (linElemSD) + CALL DEALLOCATE (coefVar) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER +END SUBROUTINE KernelAssembleIsoDiffMat + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +END MODULE KernelAssembleDiffusionMatrix_Method diff --git a/src/submodules/AbstractKernel/src/AbstractKernel_Class@AssembleTanmatMethods.F90 b/src/submodules/AbstractKernel/src/AbstractKernel_Class@AssembleTanmatMethods.F90 index eca0fa303..664415488 100644 --- a/src/submodules/AbstractKernel/src/AbstractKernel_Class@AssembleTanmatMethods.F90 +++ b/src/submodules/AbstractKernel/src/AbstractKernel_Class@AssembleTanmatMethods.F90 @@ -115,6 +115,41 @@ END IF END PROCEDURE obj_AssembleStiffnessMat +!---------------------------------------------------------------------------- +! AssembleDiffusionMat +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_AssembleDiffusionMat +CHARACTER(*), PARAMETER :: myName = "obj_AssembleDiffusionMat()" +TYPE(CPUTime_) :: TypeCPUTime + +IF (obj%showTime) CALL TypeCPUTime%SetStartTime() + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + +! NOTE: this is for isotropic one +CALL KernelAssembleDiffusionMatrix(mat=obj%diffusionMat, & + & coefficient=obj%scalarCoefficient, dom=obj%dom, & + & cellFE=obj%cellFE, linCellFE=obj%linCellFE, & + & spaceElemSD=obj%spaceElemSD, linSpaceElemSD=obj%linSpaceElemSD, & + & reset=.TRUE.) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER + +IF (obj%showTime) THEN + CALL TypeCPUTime%SetEndTime() + CALL obj%showTimeFile%WRITE(val=TypeCPUTime%GetStringForKernelLog( & + & currentTime=obj%currentTime, currentTimeStep=obj%currentTimeStep, & + & methodName=myName)) +END IF +END PROCEDURE obj_AssembleDiffusionMat + !---------------------------------------------------------------------------- ! AssembleDampingMatrix !---------------------------------------------------------------------------- From b28d444ac40d9654b56214bac62024aff0076fa7 Mon Sep 17 00:00:00 2001 From: shion <106575883+shishiousan@users.noreply.github.com> Date: Sat, 6 Jan 2024 22:45:06 +0900 Subject: [PATCH 2/5] Updates in Kernel - working in progress --- .../src/AbstractKernel_Class.F90 | 37 +++++++- .../KernelAssembleDiffusionMatrix_Method.F90 | 3 +- .../AbstractKernel/src/KernelUtility.F90 | 1 + .../AbstractKernel_Class@MaterialMethods.F90 | 85 +++++++++++++++++++ 4 files changed, 123 insertions(+), 3 deletions(-) diff --git a/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 b/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 index 621d443ea..a2cc619e0 100644 --- a/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 +++ b/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 @@ -416,7 +416,7 @@ MODULE AbstractKernel_Class CLASS(MatrixField_), POINTER :: stiffnessMat => NULL() !! Global Stiffness matrix CLASS(MatrixField_), POINTER :: diffusionMat => NULL() - !! Global Stiffness matrix + !! Global diffusion matrix CLASS(MatrixField_), POINTER :: massMat => NULL() !! Global mass matrix CLASS(MatrixField_), POINTER :: dampingMat => NULL() @@ -462,7 +462,7 @@ MODULE AbstractKernel_Class TYPE(AbstractVectorMeshFieldPointer_), ALLOCATABLE :: strain(:) !! Strain tensor !! This will be a tensor mesh field - ! TYPE(AbstractScalarMeshFieldPointer_), ALLOCATABLE :: phase_velocity(:) + ! TYPE(AbstractScalarMeshFieldPointer_), ALLOCATABLE, target:: phase_velocity(:) !! phase_velocity !! This will be a scalar mesh field TYPE(AbstractScalarMeshFieldPointer_), ALLOCATABLE :: scalarCoefficient(:) @@ -645,6 +645,12 @@ MODULE AbstractKernel_Class PROCEDURE, PUBLIC, PASS(obj) :: SetDampingProperties => & & obj_SetDampingProperties !! Set Lame parameters for isotropic elasticity + PROCEDURE, PUBLIC, PASS(obj) :: InitiateScalarCoefficient => & + & obj_InitiateScalarCoefficient + + PROCEDURE, PUBLIC, PASS(obj) :: SetScalarCoefficient => & + & obj_SetScalarCoefficient + PROCEDURE, PUBLIC, PASS(obj) :: SetMaterialToDomain => & & obj_SetMaterialToDomain !! Set material to mesh @@ -1564,6 +1570,33 @@ MODULE SUBROUTINE obj_SetDampingProperties(obj) END SUBROUTINE obj_SetDampingProperties END INTERFACE +!---------------------------------------------------------------------------- +! obj_InitiateScalarCoefficient@MaterialMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-01-06 +! summary: Initiate scalar coefficient for diffusion matrix + +INTERFACE + MODULE SUBROUTINE obj_InitiateScalarCoefficient(obj) + CLASS(AbstractKernel_), INTENT(INOUT) :: obj + END SUBROUTINE obj_InitiateScalarCoefficient +END INTERFACE + +!---------------------------------------------------------------------------- +! obj_SetScalarCoefficient@MaterialMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2024-01-06 +! summary: Set scalar coefficient for diffusion matrix + +INTERFACE + MODULE SUBROUTINE obj_SetScalarCoefficient(obj) + CLASS(AbstractKernel_), INTENT(INOUT) :: obj + END SUBROUTINE obj_SetScalarCoefficient +END INTERFACE !---------------------------------------------------------------------------- ! AddDirichletBC@BCMethods !---------------------------------------------------------------------------- diff --git a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 index c9fb5ee6f..f6924689e 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 @@ -132,7 +132,8 @@ SUBROUTINE KernelAssembleIsoDiffMat(mat, coefficient, dom, cellFE, & CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & & geoElemSD=linElemSD) - Mmat = DiffusionMatrix(test=elemsd, trial=elemsd, k=coefVar, krank=coefVar) + Mmat = DiffusionMatrix(test=elemsd, trial=elemsd, k=coefVar, & + & krank=TypeFEVariableScalar) CALL mat%Set(globalNode=nptrs, VALUE=Mmat, & & storageFMT=DOF_FMT, scale=1.0_DFP, addContribution=.TRUE.) diff --git a/src/modules/AbstractKernel/src/KernelUtility.F90 b/src/modules/AbstractKernel/src/KernelUtility.F90 index c1b5c5866..32bba7673 100644 --- a/src/modules/AbstractKernel/src/KernelUtility.F90 +++ b/src/modules/AbstractKernel/src/KernelUtility.F90 @@ -18,6 +18,7 @@ MODULE KernelUtility USE KernelAssembleMassMatrix_Method USE KernelAssembleDampingMatrix_Method USE KernelAssembleStiffnessMatrix_Method +USE KernelAssembleDiffusionMatrix_Method USE KernelAssembleElastoDynaMatrix_Method USE KernelAssembleSTElastoDynaMatrix_Method USE KernelAssembleBodySource_Method diff --git a/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 b/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 index d6af229ff..878b6995a 100644 --- a/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 +++ b/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 @@ -300,6 +300,59 @@ & 'child classes') END PROCEDURE obj_InitiateMaterialProperties +!---------------------------------------------------------------------------- +! InitiateScalarCoefficient +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_InitiateScalarCoefficient +CHARACTER(*), PARAMETER :: myName = "obj_InitiateScalarCoefficient()" +LOGICAL(LGT) :: isok +INTEGER(I4B) :: ii, tsize +TYPE(CPUTime_) :: TypeCPUTime + +IF (obj%showTime) CALL TypeCPUTime%SetStartTime() + +#ifdef DEBUG_VER +CALL e%raiseInformation(modName//'::'//myName//' - '// & + & '[START]') +#endif + +isok = ALLOCATED(obj%solidMaterial) +IF (.NOT. isok) THEN + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[NOTHING TODO] :: AbstractKernel_::obj%solidMaterial is not allocated.') + RETURN +END IF + +isok = ASSOCIATED(obj%dom) +IF (.NOT. isok) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: AbstractKernel_::obj%dom not ASSOCIATED.') + RETURN +END IF + +tsize = obj%dom%GetTotalMesh(dim=obj%nsd) +ALLOCATE (obj%scalarCoefficient(tsize)) +DO ii = 1, tsize; obj%scalarCoefficient(ii)%ptr => NULL(); END DO + +CALL KernelInitiateScalarProperty(vars=obj%scalarCoefficient, & + & materials=obj%solidMaterial, dom=obj%dom, nnt=obj%nnt, & + & varname="scalarCoefficient", matid=obj%SOLID_MATERIAL_ID, & + & engine=obj%engine%chars()) + +#ifdef DEBUG_VER +CALL e%raiseInformation(modName//'::'//myName//' - '// & + & '[END]') +#endif + +IF (obj%showTime) THEN + CALL TypeCPUTime%SetEndTime() + CALL obj%showTimeFile%WRITE(val=TypeCPUTime%GetStringForKernelLog( & + & currentTime=obj%currentTime, currentTimeStep=obj%currentTimeStep, & + & methodName=myName)) +END IF +END PROCEDURE obj_InitiateScalarCoefficient + !---------------------------------------------------------------------------- ! SetMassDensity !---------------------------------------------------------------------------- @@ -408,6 +461,38 @@ END IF END PROCEDURE obj_SetDampingProperties +!---------------------------------------------------------------------------- +! SetScalarCoefficient +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_SetScalarCoefficient +CHARACTER(*), PARAMETER :: myName = "obj_SetScalarCoefficient()" +TYPE(CPUTime_) :: TypeCPUTime + +IF (obj%showTime) CALL TypeCPUTime%SetStartTime() + +#ifdef DEBUG_VER +CALL e%raiseInformation(modName//'::'//myName//' - '// & + & '[START]') +#endif + +CALL KernelSetScalarProperty(vars=obj%scalarCoefficient, & + & materials=obj%solidMaterial, dom=obj%dom, times=obj%timeVec, & + & varname="scalarCoefficient", matid=obj%SOLID_MATERIAL_ID) + +#ifdef DEBUG_VER +CALL e%raiseInformation(modName//'::'//myName//' - '// & + & '[END]') +#endif + +IF (obj%showTime) THEN + CALL TypeCPUTime%SetEndTime() + CALL obj%showTimeFile%WRITE(val=TypeCPUTime%GetStringForKernelLog( & + & currentTime=obj%currentTime, currentTimeStep=obj%currentTimeStep, & + & methodName=myName)) +END IF +END PROCEDURE obj_SetScalarCoefficient + !---------------------------------------------------------------------------- ! SetConstantMatProps !---------------------------------------------------------------------------- From 8a2ef40b5d04313cd3e3b0c8f834d6761a6af1bb Mon Sep 17 00:00:00 2001 From: shion <106575883+shishiousan@users.noreply.github.com> Date: Tue, 9 Jan 2024 15:50:42 +0900 Subject: [PATCH 3/5] Updates in AbstractKernel_Class - working in progress --- .../src/AbstractKernel_Class.F90 | 6 +- .../src/KernelAssembleBodySource_Method.F90 | 200 +++++++++++++ .../src/KernelAssemblePointSource_Method.F90 | 73 +++++ .../KernelAssembleSurfaceSource_Method.F90 | 148 ++++++++++ .../AbstractKernel_Class@ApplyICMethods.F90 | 277 ++++++++++++++---- ...AbstractKernel_Class@ImportTomlMethods.F90 | 7 +- ...ractKernel_Class@InitiateFieldsMethods.F90 | 20 +- .../AbstractKernel_Class@MaterialMethods.F90 | 18 +- 8 files changed, 682 insertions(+), 67 deletions(-) diff --git a/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 b/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 index a2cc619e0..5f3e0ce5f 100644 --- a/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 +++ b/src/modules/AbstractKernel/src/AbstractKernel_Class.F90 @@ -1579,8 +1579,9 @@ END SUBROUTINE obj_SetDampingProperties ! summary: Initiate scalar coefficient for diffusion matrix INTERFACE - MODULE SUBROUTINE obj_InitiateScalarCoefficient(obj) + MODULE SUBROUTINE obj_InitiateScalarCoefficient(obj, varname) CLASS(AbstractKernel_), INTENT(INOUT) :: obj + CHARACTER(*), OPTIONAL, INTENT(IN) :: varname END SUBROUTINE obj_InitiateScalarCoefficient END INTERFACE @@ -1593,8 +1594,9 @@ END SUBROUTINE obj_InitiateScalarCoefficient ! summary: Set scalar coefficient for diffusion matrix INTERFACE - MODULE SUBROUTINE obj_SetScalarCoefficient(obj) + MODULE SUBROUTINE obj_SetScalarCoefficient(obj, varname) CLASS(AbstractKernel_), INTENT(INOUT) :: obj + CHARACTER(*), OPTIONAL, INTENT(IN) :: varname END SUBROUTINE obj_SetScalarCoefficient END INTERFACE !---------------------------------------------------------------------------- diff --git a/src/modules/AbstractKernel/src/KernelAssembleBodySource_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleBodySource_Method.F90 index fc4173309..14ad0d6da 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleBodySource_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleBodySource_Method.F90 @@ -37,6 +37,8 @@ MODULE KernelAssembleBodySource_Method INTERFACE KernelAssembleBodySource MODULE PROCEDURE KernelAssembleBodySource1 MODULE PROCEDURE KernelAssembleBodySource2 + MODULE PROCEDURE KernelAssembleBodySource3 + MODULE PROCEDURE KernelAssembleBodySource4 END INTERFACE KernelAssembleBodySource CONTAINS @@ -191,6 +193,7 @@ SUBROUTINE KernelAssembleBodySource2(rhs, dom, bodyVec, cellFE, & CALL Reallocate(xij, nsd, nns) CALL Reallocate(fevec, nsd, nns) + CALL Reallocate(forceVec, nsd, nns) DO iel = meshptr%minElemNum, meshptr%maxElemNum @@ -233,4 +236,201 @@ SUBROUTINE KernelAssembleBodySource2(rhs, dom, bodyVec, cellFE, & END SUBROUTINE KernelAssembleBodySource2 +!---------------------------------------------------------------------------- +! KernelAssembleBodySource +!---------------------------------------------------------------------------- + +SUBROUTINE KernelAssembleBodySource3(rhs, dom, func, cellFE, & + & linCellFE, spaceElemSD, linSpaceElemSD, reset, scale, times) + CLASS(ScalarField_), INTENT(INOUT) :: rhs + CLASS(Domain_), INTENT(INOUT) :: dom + CLASS(UserFunction_), INTENT(INOUT) :: func + TYPE(FiniteElementPointer_), INTENT(INOUT) :: cellFE(:) + TYPE(FiniteElementPointer_), INTENT(INOUT) :: linCellFE(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: spaceElemSD(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: linSpaceElemSD(:) + LOGICAL(LGT), INTENT(IN) :: reset + REAL(DFP), INTENT(IN) :: scale + REAL(DFP), OPTIONAL, INTENT(IN) :: times(:) + + ! internal variables + CHARACTER(*), PARAMETER :: myName = "KernelAssembleBodySource3()" + LOGICAL(LGT) :: problem, isok + INTEGER(I4B) :: tmesh, nsd, id, nns, iel + INTEGER(I4B), ALLOCATABLE :: nptrs(:) + REAL(DFP), ALLOCATABLE :: fevec(:), xij(:, :) + TYPE(FEVariable_) :: forceVar + TYPE(ElemShapeData_) :: elemsd, linElemSD + CLASS(Mesh_), POINTER :: meshptr + CLASS(ReferenceElement_), POINTER :: refelem + CLASS(FiniteElement_), POINTER :: spaceFE, linSpaceFE + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + + IF (reset) CALL rhs%Set(VALUE=0.0_DFP) + + nsd = dom%GetNSD() + tmesh = dom%GetTotalMesh(dim=nsd) + NULLIFY (meshptr, refelem, spaceFE, linSpaceFE) + + DO id = 1, tmesh + meshptr => dom%GetMeshPointer(dim=nsd, entityNum=id) + problem = meshptr%isEmpty() + IF (problem) CYCLE + + spaceFE => cellFE(id)%ptr + linSpaceFE => linCellFE(id)%ptr + + elemsd = spaceElemSD(id) + linElemSD = linSpaceElemSD(id) + + refelem => meshptr%GetRefElemPointer() + nns = (.nne.refelem) + + CALL Reallocate(xij, nsd, nns) + CALL Reallocate(fevec, nns) + + DO iel = meshptr%minElemNum, meshptr%maxElemNum + + isok = meshptr%isElementPresent(iel) + IF (.NOT. isok) CYCLE + + nptrs = meshptr%GetConnectivity(iel) + + CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) + + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + & geoElemSD=linElemSD) + + CALL func%Get(fevar=forceVar, xij=xij, times=times) + + fevec = ForceVector(test=elemsd, c=forceVar, & + & crank=TypeFEVariableScalar) + + CALL rhs%Set(globalNode=nptrs, scale=scale, & + & addContribution=.TRUE., VALUE=fevec) + + END DO + + END DO + + NULLIFY (meshptr, refelem, spaceFE, linSpaceFE) + IF (ALLOCATED(nptrs)) DEALLOCATE (nptrs) + IF (ALLOCATED(fevec)) DEALLOCATE (fevec) + CALL DEALLOCATE (forceVar) + CALL DEALLOCATE (elemsd) + CALL DEALLOCATE (linElemSD) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER + +END SUBROUTINE KernelAssembleBodySource3 + +!---------------------------------------------------------------------------- +! KernelAssembleBodySource +!---------------------------------------------------------------------------- + +SUBROUTINE KernelAssembleBodySource4(rhs, dom, bodyVec, cellFE, & + & linCellFE, spaceElemSD, linSpaceElemSD, reset, scale) + CLASS(ScalarField_), INTENT(INOUT) :: rhs + CLASS(Domain_), INTENT(INOUT) :: dom + CLASS(ScalarField_), INTENT(INOUT) :: bodyVec + TYPE(FiniteElementPointer_), INTENT(INOUT) :: cellFE(:) + TYPE(FiniteElementPointer_), INTENT(INOUT) :: linCellFE(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: spaceElemSD(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: linSpaceElemSD(:) + LOGICAL(LGT), INTENT(IN) :: reset + REAL(DFP), INTENT(IN) :: scale + + ! internal variables + CHARACTER(*), PARAMETER :: myName = "KernelAssembleBodySource4()" + LOGICAL(LGT) :: problem, isok + INTEGER(I4B) :: tmesh, nsd, id, nns, iel + INTEGER(I4B), ALLOCATABLE :: nptrs(:) + REAL(DFP), ALLOCATABLE :: fevec(:), xij(:, :), forceVec(:) + TYPE(FEVariable_) :: forceVar + TYPE(ElemShapeData_) :: elemsd, linElemSD + CLASS(Mesh_), POINTER :: meshptr + CLASS(ReferenceElement_), POINTER :: refelem + CLASS(FiniteElement_), POINTER :: spaceFE, linSpaceFE + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + + IF (reset) CALL rhs%Set(VALUE=0.0_DFP) + + nsd = dom%GetNSD() + tmesh = dom%GetTotalMesh(dim=nsd) + NULLIFY (meshptr, refelem, spaceFE, linSpaceFE) + + DO id = 1, tmesh + meshptr => dom%GetMeshPointer(dim=nsd, entityNum=id) + problem = meshptr%isEmpty() + IF (problem) CYCLE + + spaceFE => cellFE(id)%ptr + linSpaceFE => linCellFE(id)%ptr + + elemsd = spaceElemSD(id) + linElemSD = linSpaceElemSD(id) + + refelem => meshptr%GetRefElemPointer() + nns = (.nne.refelem) + + CALL Reallocate(xij, nsd, nns) + CALL Reallocate(fevec, nns) + CALL Reallocate(forceVec, nns) + + DO iel = meshptr%minElemNum, meshptr%maxElemNum + + isok = meshptr%isElementPresent(iel) + IF (.NOT. isok) CYCLE + + nptrs = meshptr%GetConnectivity(iel) + + CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) + + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + & geoElemSD=linElemSD) + + CALL bodyVec%Get(VALUE=forceVec, globalNode=nptrs) + + forceVar = NodalVariable(val=forceVec, rank=TypeFEVariableScalar, & + & vartype=TypeFEVariableSpace) + + fevec = ForceVector(test=elemsd, c=forceVar, & + & crank=TypeFEVariableScalar) + + CALL rhs%Set(globalNode=nptrs, scale=scale, & + & addContribution=.TRUE., VALUE=fevec) + + END DO + + END DO + + NULLIFY (meshptr, refelem, spaceFE, linSpaceFE) + IF (ALLOCATED(nptrs)) DEALLOCATE (nptrs) + IF (ALLOCATED(fevec)) DEALLOCATE (fevec) + IF (ALLOCATED(forceVec)) DEALLOCATE (forceVec) + CALL DEALLOCATE (forceVar) + CALL DEALLOCATE (elemsd) + CALL DEALLOCATE (linElemSD) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER + +END SUBROUTINE KernelAssembleBodySource4 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- END MODULE KernelAssembleBodySource_Method diff --git a/src/modules/AbstractKernel/src/KernelAssemblePointSource_Method.F90 b/src/modules/AbstractKernel/src/KernelAssemblePointSource_Method.F90 index 6794703d4..e02f976da 100644 --- a/src/modules/AbstractKernel/src/KernelAssemblePointSource_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssemblePointSource_Method.F90 @@ -36,6 +36,7 @@ MODULE KernelAssemblePointSource_Method INTERFACE KernelAssemblePointSource MODULE PROCEDURE KernelAssemblePointSource1 + MODULE PROCEDURE KernelAssemblePointSource2 END INTERFACE KernelAssemblePointSource CONTAINS @@ -113,4 +114,76 @@ SUBROUTINE KernelAssemblePointSource1(rhs, nbcPtrs, reset, scale, times) END SUBROUTINE KernelAssemblePointSource1 +!---------------------------------------------------------------------------- +! KernelAssemblePointSource +!---------------------------------------------------------------------------- + +SUBROUTINE KernelAssemblePointSource2(rhs, nbcPtrs, reset, scale, times) + CLASS(ScalarField_), INTENT(INOUT) :: rhs + TYPE(NeumannBCPointer_), INTENT(INOUT) :: nbcPtrs(:) + LOGICAL(LGT), INTENT(IN) :: reset + REAL(DFP), INTENT(IN) :: scale + REAL(DFP), OPTIONAL, INTENT(IN) :: times(:) + + ! internal variables + CHARACTER(*), PARAMETER :: myName = "KernelAssemblePointSource2()" + CLASS(NeumannBC_), POINTER :: nbc + LOGICAL(LGT) :: problem, isSelectionByNodeNum, istimes + INTEGER(I4B) :: tnbc, nbcNo, idof + INTEGER(I4B), ALLOCATABLE :: nodeNum(:) + REAL(DFP), ALLOCATABLE :: nodalValue(:, :) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + + IF (reset) CALL rhs%Set(VALUE=0.0_DFP) + + istimes = PRESENT(times) + IF (istimes) THEN + nbcNo = SIZE(times) + problem = nbcNo .NE. 1_I4B + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: size of times should be 1.') + RETURN + END IF + END IF + + tnbc = SIZE(nbcPtrs) + NULLIFY (nbc) + + DO nbcNo = 1, tnbc + nbc => NULL() + nbc => nbcPtrs(nbcNo)%ptr + + problem = .NOT. ASSOCIATED(nbc) + IF (problem) CYCLE + CALL nbc%GetParam(isSelectionByNodeNum=isSelectionByNodeNum) + problem = .NOT. isSelectionByNodeNum + IF (problem) CYCLE + + CALL nbc%Get(nodeNum=nodeNum, nodalValue=nodalValue, times=times) + idof = nbc%GetDOFNo() + + CALL rhs%Set(globalNode=nodeNum, VALUE=nodalValue(:, 1), & + & scale=scale, addContribution=.TRUE.) + END DO + + NULLIFY (nbc) + IF (ALLOCATED(nodeNum)) DEALLOCATE (nodeNum) + IF (ALLOCATED(nodalValue)) DEALLOCATE (nodalValue) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER + +END SUBROUTINE KernelAssemblePointSource2 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + END MODULE KernelAssemblePointSource_Method diff --git a/src/modules/AbstractKernel/src/KernelAssembleSurfaceSource_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleSurfaceSource_Method.F90 index e4e1a3617..6ac36ce79 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleSurfaceSource_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleSurfaceSource_Method.F90 @@ -39,6 +39,7 @@ MODULE KernelAssembleSurfaceSource_Method MODULE PROCEDURE KernelAssembleSurfaceSource1 MODULE PROCEDURE KernelAssembleSurfaceSource2 MODULE PROCEDURE KernelAssembleSurfaceSource3 + MODULE PROCEDURE KernelAssembleSurfaceSource4 END INTERFACE KernelAssembleSurfaceSource CONTAINS @@ -391,4 +392,151 @@ SUBROUTINE KernelAssembleSurfaceSource3(rhs, extField, dom, nbcPtrs, fe, & END SUBROUTINE KernelAssembleSurfaceSource3 +!---------------------------------------------------------------------------- +! KernelAssembleSurfaceSource +!---------------------------------------------------------------------------- + +SUBROUTINE KernelAssembleSurfaceSource4(rhs, extField, dom, nbcPtrs, fe, & + & linFE, spaceElemSD, linSpaceElemSD, reset, scale, times) + CLASS(ScalarField_), INTENT(INOUT) :: rhs + CLASS(ScalarField_), INTENT(INOUT) :: extField + CLASS(Domain_), INTENT(INOUT) :: dom + TYPE(NeumannBCPointer_), INTENT(INOUT) :: nbcPtrs(:) + TYPE(FiniteElementPointer_), INTENT(INOUT) :: fe(:) + TYPE(FiniteElementPointer_), INTENT(INOUT) :: linFE(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: spaceElemSD(:) + TYPE(ElemShapeData_), INTENT(INOUT) :: linSpaceElemSD(:) + LOGICAL(LGT), INTENT(IN) :: reset + REAL(DFP), INTENT(IN) :: scale + REAL(DFP), OPTIONAL, INTENT(IN) :: times(:) + + ! internal variables + CHARACTER(*), PARAMETER :: myName = "KernelAssembleSurfaceSource4()" + TYPE(FEVariable_) :: forceVar + TYPE(ElemShapeData_) :: elemsd, linElemSD + CLASS(Mesh_), POINTER :: meshptr + CLASS(ReferenceElement_), POINTER :: refelem + CLASS(FiniteElement_), POINTER :: spaceFE, linSpaceFE + CLASS(NeumannBC_), POINTER :: nbc + LOGICAL(LGT) :: problem, isNormal, isTangent, isSelectionByMeshID + INTEGER(I4B) :: tmesh, nsd, id, nns, iel, tnbc, nbcNo, idof, jd + INTEGER(I4B), ALLOCATABLE :: nptrs(:), meshID(:) + REAL(DFP), ALLOCATABLE :: fevec(:), xij(:, :), forceVec(:) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + + IF (reset) CALL rhs%Set(VALUE=0.0_DFP) + + tnbc = SIZE(nbcPtrs) + nsd = dom%GetNSD() + NULLIFY (meshptr, refelem, spaceFE, linSpaceFE, nbc) + + DO nbcNo = 1, tnbc + nbc => nbcPtrs(nbcNo)%ptr + problem = .NOT. ASSOCIATED(nbc) + IF (problem) CYCLE + + CALL nbc%GetParam(useExternal=problem) + IF (problem) CYCLE + + CALL extField%ApplyDirichletBC(dbc=nbc, times=times) + nbc => NULL() + END DO + + DO nbcNo = 1, tnbc + nbc => nbcPtrs(nbcNo)%ptr + problem = .NOT. ASSOCIATED(nbc) + IF (problem) CYCLE + + CALL nbc%GetParam(isSelectionByMeshID=isSelectionByMeshID, & + & isNormal=isNormal, isTangent=isTangent) + + problem = .NOT. isSelectionByMeshID + IF (problem) THEN + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[SKIPPING] :: Currently, found isSelectionByMeshID false.') + CYCLE + END IF + + problem = isTangent .OR. isNormal + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: Currently, normal and tangential '// & + & ' Neumann boundary condition are not supported.') + RETURN + END IF + + meshID = nbc%GetMeshID(dim=nsd - 1_I4B) + idof = nbc%GetDOFNo() + tmesh = SIZE(meshID) + + DO id = 1, tmesh + jd = meshID(id) + meshptr => dom%GetMeshPointer(dim=nsd - 1_I4B, entityNum=jd) + + problem = meshptr%isEmpty() + IF (problem) CYCLE + + spaceFE => fe(jd)%ptr + linSpaceFE => linFE(jd)%ptr + + elemsd = spaceElemSD(jd) + linElemSD = linSpaceElemSD(jd) + + refelem => meshptr%GetRefElemPointer() + nns = (.nne.refelem) + + CALL Reallocate(xij, nsd, nns) + CALL Reallocate(fevec, nns) + CALL Reallocate(forceVec, nns) + + DO iel = meshptr%minElemNum, meshptr%maxElemNum + + problem = .NOT. meshptr%IsElementPresent(iel) + IF (problem) CYCLE + + nptrs = meshptr%GetConnectivity(iel) + CALL dom%GetNodeCoord(nodeCoord=xij, globalNode=nptrs) + + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + & geoElemSD=linElemSD) + + CALL extField%Get(VALUE=forceVec, globalNode=nptrs) + + forceVar = NodalVariable(val=forceVec, rank=TypeFEVariableScalar, & + & vartype=TypeFEVariableSpace) + + fevec = ForceVector(test=elemsd, c=forceVar, & + & crank=TypeFEVariableScalar) + + CALL rhs%Set(globalNode=nptrs, VALUE=fevec, scale=scale, & + & addContribution=.TRUE.) + END DO + END DO + + END DO + + NULLIFY (meshptr, refelem, spaceFE, linSpaceFE, nbc) + + IF (ALLOCATED(nptrs)) DEALLOCATE (nptrs) + IF (ALLOCATED(meshID)) DEALLOCATE (meshID) + IF (ALLOCATED(fevec)) DEALLOCATE (fevec) + IF (ALLOCATED(xij)) DEALLOCATE (xij) + + CALL DEALLOCATE (forceVar) + CALL DEALLOCATE (elemsd) + CALL DEALLOCATE (linElemSD) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER + +END SUBROUTINE KernelAssembleSurfaceSource4 +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- END MODULE KernelAssembleSurfaceSource_Method diff --git a/src/submodules/AbstractKernel/src/AbstractKernel_Class@ApplyICMethods.F90 b/src/submodules/AbstractKernel/src/AbstractKernel_Class@ApplyICMethods.F90 index e0e3351c0..7b031c5aa 100644 --- a/src/submodules/AbstractKernel/src/AbstractKernel_Class@ApplyICMethods.F90 +++ b/src/submodules/AbstractKernel/src/AbstractKernel_Class@ApplyICMethods.F90 @@ -21,93 +21,251 @@ CONTAINS !---------------------------------------------------------------------------- -! ApplyIC +! !---------------------------------------------------------------------------- -MODULE PROCEDURE obj_ApplyIC -CHARACTER(*), PARAMETER :: myName = "obj_ApplyIC()" -CHARACTER(:), ALLOCATABLE :: name0 -LOGICAL(LGT) :: problem, isfunc, isext -TYPE(CPUTime_) :: TypeCPUTime - +SUBROUTINE applyIC2Vec(obj, name, func, extField, times, ivar, & + & idof, spaceCompo, timeCompo) + CLASS(AbstractKernel_), INTENT(INOUT) :: obj + CHARACTER(*), OPTIONAL, INTENT(IN) :: name + CLASS(UserFunction_), OPTIONAL, INTENT(INOUT) :: func + CLASS(AbstractNodeField_), OPTIONAL, INTENT(INOUT) :: extField + REAL(DFP), OPTIONAL, INTENT(IN) :: times(:) + INTEGER(I4B), OPTIONAL, INTENT(IN) :: ivar + INTEGER(I4B), OPTIONAL, INTENT(IN) :: idof + INTEGER(I4B), OPTIONAL, INTENT(IN) :: spaceCompo + INTEGER(I4B), OPTIONAL, INTENT(IN) :: timeCompo -IF (obj%showTime) CALL TypeCPUTime%SetStartTime() + CHARACTER(*), PARAMETER :: myName = "ApplyIC2Vec()" + CHARACTER(:), ALLOCATABLE :: name0 + LOGICAL(LGT) :: problem, isfunc, isext #ifdef DEBUG_VER -CALL e%RaiseInformation(modName//'::'//myName//' - '// & - & '[START] ') + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') #endif DEBUG_VER -IF (PRESENT(name)) THEN - name0 = UpperCase(name) -ELSE - name0 = "NONE" -END IF + IF (PRESENT(name)) THEN + name0 = UpperCase(name) + ELSE + name0 = "NONE" + END IF -isext = PRESENT(extField) -isfunc = PRESENT(func) + isext = PRESENT(extField) + isfunc = PRESENT(func) -SELECT CASE (name0) -CASE ("DISP", "DISPLACEMENT") - problem = .NOT. ASSOCIATED(obj%displacement) - IF (problem) THEN - CALL e%RaiseError(modName//'::'//myName//' - '// & - & '[INTERNAL ERROR] :: displacement is not ASSOCIATED.') - RETURN - END IF + SELECT CASE (name0) + CASE ("DISP", "DISPLACEMENT") + problem = .NOT. ASSOCIATED(obj%displacement) + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: displacement is not ASSOCIATED.') + RETURN + END IF - CALL obj%displacement%Set(VALUE=0.0_DFP) + CALL obj%displacement%Set(VALUE=0.0_DFP) - IF (isfunc) THEN - CALL obj%displacement%Set(func=func, times=times, & - & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) - END IF + IF (isfunc) THEN + CALL obj%displacement%Set(func=func, times=times, & + & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + END IF - IF (isext) THEN - CALL obj%displacement%Copy(obj2=extField) - END IF + IF (isext) THEN + CALL obj%displacement%Copy(obj2=extField) + END IF + + CASE ("VEL", "VELOCITY") + problem = .NOT. ASSOCIATED(obj%velocity) + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: velocity is not ASSOCIATED.') + RETURN + END IF + + CALL obj%velocity%Set(VALUE=0.0_DFP) + + IF (isfunc) THEN + CALL obj%velocity%Set(func=func, times=times, & + & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + END IF -CASE ("VEL", "VELOCITY") - problem = .NOT. ASSOCIATED(obj%velocity) - IF (problem) THEN + IF (isext) THEN + CALL obj%velocity%Copy(obj2=extField) + END IF + + CASE ("ACC", "ACCELERATION") + problem = .NOT. ASSOCIATED(obj%acceleration) + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: acceleration is not ASSOCIATED.') + RETURN + END IF + + CALL obj%acceleration%Set(VALUE=0.0_DFP) + + IF (isfunc) THEN + CALL obj%acceleration%Set(func=func, times=times, & + & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + END IF + + IF (isext) THEN + CALL obj%acceleration%Copy(obj2=extField) + END IF + + CASE DEFAULT CALL e%RaiseError(modName//'::'//myName//' - '// & - & '[INTERNAL ERROR] :: velocity is not ASSOCIATED.') + & '[INTERNAL ERROR] :: No case found for name0') RETURN - END IF + END SELECT - CALL obj%velocity%Set(VALUE=0.0_DFP) +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER - IF (isfunc) THEN - CALL obj%velocity%Set(func=func, times=times, & - & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) - END IF +END SUBROUTINE applyIC2Vec - IF (isext) THEN - CALL obj%velocity%Copy(obj2=extField) +!---------------------------------------------------------------------------- +! ApplyIC2Scalar +!---------------------------------------------------------------------------- + +SUBROUTINE applyIC2Scalar(obj, name, func, extField, times, ivar, & + & idof, spaceCompo, timeCompo) + CLASS(AbstractKernel_), INTENT(INOUT) :: obj + CHARACTER(*), OPTIONAL, INTENT(IN) :: name + CLASS(UserFunction_), OPTIONAL, INTENT(INOUT) :: func + CLASS(AbstractNodeField_), OPTIONAL, INTENT(INOUT) :: extField + REAL(DFP), OPTIONAL, INTENT(IN) :: times(:) + INTEGER(I4B), OPTIONAL, INTENT(IN) :: ivar + INTEGER(I4B), OPTIONAL, INTENT(IN) :: idof + INTEGER(I4B), OPTIONAL, INTENT(IN) :: spaceCompo + INTEGER(I4B), OPTIONAL, INTENT(IN) :: timeCompo + + CHARACTER(*), PARAMETER :: myName = "applyIC2Scalar()" + CHARACTER(:), ALLOCATABLE :: name0 + LOGICAL(LGT) :: problem, isfunc, isext + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + + IF (PRESENT(name)) THEN + name0 = UpperCase(name) + ELSE + name0 = "NONE" END IF -CASE ("ACC", "ACCELERATION") - problem = .NOT. ASSOCIATED(obj%acceleration) - IF (problem) THEN + isext = PRESENT(extField) + isfunc = PRESENT(func) + + SELECT CASE (name0) + CASE ("PRESSURE", "PRES") + problem = .NOT. ASSOCIATED(obj%pressure) + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: pressure is not ASSOCIATED.') + RETURN + END IF + + CALL obj%pressure%Set(VALUE=0.0_DFP) + + IF (isfunc) THEN + CALL obj%pressure%Set(func=func, times=times, & + & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + END IF + + IF (isext) THEN + CALL obj%pressure%Copy(obj2=extField) + END IF + + CASE ("PVEL", "PVELOCITY") + problem = .NOT. ASSOCIATED(obj%p_velocity) + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: p_velocity is not ASSOCIATED.') + RETURN + END IF + + CALL obj%p_velocity%Set(VALUE=0.0_DFP) + + IF (isfunc) THEN + CALL obj%p_velocity%Set(func=func, times=times, & + & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + END IF + + IF (isext) THEN + CALL obj%p_velocity%Copy(obj2=extField) + END IF + + CASE ("PACC", "PACCELERATION") + problem = .NOT. ASSOCIATED(obj%p_acceleration) + IF (problem) THEN + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: p_acceleration is not ASSOCIATED.') + RETURN + END IF + + CALL obj%p_acceleration%Set(VALUE=0.0_DFP) + + IF (isfunc) THEN + CALL obj%p_acceleration%Set(func=func, times=times, & + & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + END IF + + IF (isext) THEN + CALL obj%p_acceleration%Copy(obj2=extField) + END IF + + CASE DEFAULT CALL e%RaiseError(modName//'::'//myName//' - '// & - & '[INTERNAL ERROR] :: acceleration is not ASSOCIATED.') + & '[INTERNAL ERROR] :: No case found for name0') RETURN - END IF + END SELECT - CALL obj%acceleration%Set(VALUE=0.0_DFP) +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[END] ') +#endif DEBUG_VER - IF (isfunc) THEN - CALL obj%acceleration%Set(func=func, times=times, & - & ivar=ivar, idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) - END IF +END SUBROUTINE applyIC2Scalar - IF (isext) THEN - CALL obj%acceleration%Copy(obj2=extField) - END IF +!---------------------------------------------------------------------------- +! ApplyIC +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplyIC +CHARACTER(*), PARAMETER :: myName = "obj_ApplyIC()" +TYPE(CPUTime_) :: TypeCPUTime + +IF (obj%showTime) CALL TypeCPUTime%SetStartTime() + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + & '[START] ') +#endif DEBUG_VER + +SELECT CASE (obj%problemType) +CASE (KernelProblemType%scalar) + CALL ApplyIC2Scalar(obj=obj, name=name, func=func, & + & extField=extField, times=times, ivar=ivar, & + & idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + +CASE (KernelProblemType%vector) + + CALL ApplyIC2Vec(obj=obj, name=name, func=func, & + & extField=extField, times=times, ivar=ivar, & + & idof=idof, spaceCompo=spaceCompo, timeCompo=timeCompo) + +CASE (KernelProblemType%multiPhysics) + + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[WIP] :: not implemented yet.') + RETURN CASE DEFAULT CALL e%RaiseError(modName//'::'//myName//' - '// & - & '[INTERNAL ERROR] :: No case found for name0') + & '[INTERNAL ERROR] :: No case found for KernelProblemType') RETURN END SELECT @@ -122,6 +280,7 @@ & currentTime=obj%currentTime, currentTimeStep=obj%currentTimeStep, & & methodName=myName)) END IF + END PROCEDURE obj_ApplyIC END SUBMODULE ApplyICMethods diff --git a/src/submodules/AbstractKernel/src/AbstractKernel_Class@ImportTomlMethods.F90 b/src/submodules/AbstractKernel/src/AbstractKernel_Class@ImportTomlMethods.F90 index f6e285d9b..e707512a7 100644 --- a/src/submodules/AbstractKernel/src/AbstractKernel_Class@ImportTomlMethods.F90 +++ b/src/submodules/AbstractKernel/src/AbstractKernel_Class@ImportTomlMethods.F90 @@ -39,7 +39,7 @@ & baseContinuityForSpace, quadratureTypeForSpace, ipTypeForSpace, & & basisTypeForSpace, baseInterpolationForTime, baseContinuityForTime, & & quadratureTypeForTime, ipTypeForTime, basisTypeForTime, & - & problemType, tanmatProp, astr + & problemType, tanmatProp, astr, outputPath INTEGER(I4B) :: algorithm, tSolidMaterials, tDirichletBC, tWeakDirichletBC, & & tNeumannBC, tMaterialInterfaces, origin, stat, maxIter, nsd, nnt, tdof, & @@ -105,6 +105,10 @@ & tanmatProp%raw, DEFAULT_TANMAT_PROP, & & origin=origin, stat=stat) +CALL toml_get(table, "outputPath", & + & outputPath%raw, DEFAULT_OUTPUT_PATH, & + & origin=origin, stat=stat) + CALL toml_get(table, "problemType", problemType%raw, & & DEFAULT_PROBLEM_TYPE_CHAR, origin=origin, stat=stat) @@ -352,6 +356,7 @@ & rtoleranceForResidual=rtoleranceForResidual, & & atoleranceForResidual=atoleranceForResidual, & & tanmatProp=tanmatProp%chars(), & + & outputPath=outputPath%chars(), & & tOverlappedMaterials=tOverlappedMaterials, & & showTime=showTime) diff --git a/src/submodules/AbstractKernel/src/AbstractKernel_Class@InitiateFieldsMethods.F90 b/src/submodules/AbstractKernel/src/AbstractKernel_Class@InitiateFieldsMethods.F90 index 35086e597..cb17a7c43 100644 --- a/src/submodules/AbstractKernel/src/AbstractKernel_Class@InitiateFieldsMethods.F90 +++ b/src/submodules/AbstractKernel/src/AbstractKernel_Class@InitiateFieldsMethods.F90 @@ -30,6 +30,7 @@ MODULE PROCEDURE obj_InitiateTangentMatrix CHARACTER(*), PARAMETER :: myName = "obj_InitiateTangentMatrix()" LOGICAL(LGT) :: isok +INTEGER(I4B) :: nsd0 TYPE(CPUTime_) :: TypeCPUTime IF (obj%showTime) CALL TypeCPUTime%SetStartTime() @@ -48,10 +49,23 @@ & name="MATRIX") END IF +SELECT CASE (obj%problemType) +CASE (KernelProblemType%scalar) + nsd0 = 1_I4B +CASE (KernelProblemType%vector) + nsd0 = obj%nsd +CASE (KernelProblemType%multiPhysics) + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[WIP] :: not implemented yet') +CASE default + CALL e%RaiseError(modName//'::'//myName//' - '// & + & '[INTERNAL ERROR] :: no case found for KernelProblemType') +END SELECT + CALL KernelInitiateTangentMatrix(mat=obj%tanmat, & - & linsol=obj%linsol, dom=obj%dom, nsd=obj%nsd, nnt=obj%nnt, & - & engine=obj%engine%chars(), name="tanmat", & - & matrixProp=obj%tanmatProp%chars()) + & linsol=obj%linsol, dom=obj%dom, nsd=nsd0, nnt=obj%nnt, & + & engine=obj%engine%chars(), name="tanmat", & + & matrixProp=obj%tanmatProp%chars()) #ifdef DEBUG_VER CALL e%RaiseInformation(modName//'::'//myName//' - '// & diff --git a/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 b/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 index 878b6995a..ea5598fe6 100644 --- a/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 +++ b/src/submodules/AbstractKernel/src/AbstractKernel_Class@MaterialMethods.F90 @@ -308,6 +308,7 @@ CHARACTER(*), PARAMETER :: myName = "obj_InitiateScalarCoefficient()" LOGICAL(LGT) :: isok INTEGER(I4B) :: ii, tsize +CHARACTER(:), ALLOCATABLE :: varname0 TYPE(CPUTime_) :: TypeCPUTime IF (obj%showTime) CALL TypeCPUTime%SetStartTime() @@ -331,13 +332,19 @@ RETURN END IF +IF (PRESENT(varname)) THEN + varname0 = varname +ELSE + varname0 = "scalarCoefficient" +END IF + tsize = obj%dom%GetTotalMesh(dim=obj%nsd) ALLOCATE (obj%scalarCoefficient(tsize)) DO ii = 1, tsize; obj%scalarCoefficient(ii)%ptr => NULL(); END DO CALL KernelInitiateScalarProperty(vars=obj%scalarCoefficient, & & materials=obj%solidMaterial, dom=obj%dom, nnt=obj%nnt, & - & varname="scalarCoefficient", matid=obj%SOLID_MATERIAL_ID, & + & varname=varname0, matid=obj%SOLID_MATERIAL_ID, & & engine=obj%engine%chars()) #ifdef DEBUG_VER @@ -468,6 +475,7 @@ MODULE PROCEDURE obj_SetScalarCoefficient CHARACTER(*), PARAMETER :: myName = "obj_SetScalarCoefficient()" TYPE(CPUTime_) :: TypeCPUTime +CHARACTER(:), ALLOCATABLE :: varname0 IF (obj%showTime) CALL TypeCPUTime%SetStartTime() @@ -476,9 +484,15 @@ & '[START]') #endif +IF (PRESENT(varname)) THEN + varname0 = varname +ELSE + varname0 = "scalarCoefficient" +END IF + CALL KernelSetScalarProperty(vars=obj%scalarCoefficient, & & materials=obj%solidMaterial, dom=obj%dom, times=obj%timeVec, & - & varname="scalarCoefficient", matid=obj%SOLID_MATERIAL_ID) + & varname=varname0, matid=obj%SOLID_MATERIAL_ID) #ifdef DEBUG_VER CALL e%raiseInformation(modName//'::'//myName//' - '// & From 1a44eee521938c441d32d48310d6eec890045b3f Mon Sep 17 00:00:00 2001 From: shion <106575883+shishiousan@users.noreply.github.com> Date: Wed, 10 Jan 2024 16:23:21 +0900 Subject: [PATCH 4/5] Updates in AbstractKernel - working in progress for Acoustic kernel --- .../src/KernelAssembleDiffusionMatrix_Method.F90 | 5 +++-- .../AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 index f6924689e..2acceced9 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 @@ -129,8 +129,9 @@ SUBROUTINE KernelAssembleIsoDiffMat(mat, coefficient, dom, cellFE, & CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) - CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & - & geoElemSD=linElemSD) + ! CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + ! & geoElemSD=linElemSD) + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij) Mmat = DiffusionMatrix(test=elemsd, trial=elemsd, k=coefVar, & & krank=TypeFEVariableScalar) diff --git a/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 index ed7574169..b437f7b48 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 @@ -115,8 +115,9 @@ SUBROUTINE KernelAssembleMassMatrix1(mat, massDensity, dom, cellFE, & CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) - CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & - & geoElemSD=linElemSD) + ! CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + ! & geoElemSD=linElemSD) + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij) Mmat = MassMatrix(test=elemsd, trial=elemsd, opt=tdof, & & rho=fevar, rhorank=TypeFEVariableScalar) From b2178428d647106c01c802f21c5ad8d30545b5bd Mon Sep 17 00:00:00 2001 From: shion <106575883+shishiousan@users.noreply.github.com> Date: Thu, 11 Jan 2024 19:36:41 +0900 Subject: [PATCH 5/5] Updates in AbstractKernel - Minor changes --- .../src/KernelAssembleDiffusionMatrix_Method.F90 | 5 ++--- .../AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 | 5 ++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 index 2acceced9..f6924689e 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleDiffusionMatrix_Method.F90 @@ -129,9 +129,8 @@ SUBROUTINE KernelAssembleIsoDiffMat(mat, coefficient, dom, cellFE, & CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) - ! CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & - ! & geoElemSD=linElemSD) - CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij) + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + & geoElemSD=linElemSD) Mmat = DiffusionMatrix(test=elemsd, trial=elemsd, k=coefVar, & & krank=TypeFEVariableScalar) diff --git a/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 b/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 index b437f7b48..ed7574169 100644 --- a/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 +++ b/src/modules/AbstractKernel/src/KernelAssembleMassMatrix_Method.F90 @@ -115,9 +115,8 @@ SUBROUTINE KernelAssembleMassMatrix1(mat, massDensity, dom, cellFE, & CALL dom%GetNodeCoord(globalNode=nptrs, nodeCoord=xij) - ! CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & - ! & geoElemSD=linElemSD) - CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij) + CALL spaceFE%GetGlobalElemShapeData(elemsd=elemsd, xij=xij, & + & geoElemSD=linElemSD) Mmat = MassMatrix(test=elemsd, trial=elemsd, opt=tdof, & & rho=fevar, rhorank=TypeFEVariableScalar)