Merge pull request #371 from easifem-fortran/acoustic

Updates in AbstractKernel
easifem · Jan 15, 2024 · 55efdf3 · 55efdf3
2 parents 8b44f09 + b217842
commit 55efdf3
Show file tree

Hide file tree

Showing 12 changed files with 1,026 additions and 63 deletions.
diff --git 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
@@ -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 diffusion 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, target:: 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
 
@@ -631,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
@@ -678,6 +698,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 =>  &
@@ -1548,6 +1570,35 @@ 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, varname)
+    CLASS(AbstractKernel_), INTENT(INOUT) :: obj
+    CHARACTER(*), OPTIONAL, INTENT(IN) :: varname
+  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, varname)
+    CLASS(AbstractKernel_), INTENT(INOUT) :: obj
+    CHARACTER(*), OPTIONAL, INTENT(IN) :: varname
+  END SUBROUTINE obj_SetScalarCoefficient
+END INTERFACE
 !----------------------------------------------------------------------------
 !                                                   AddDirichletBC@BCMethods
 !----------------------------------------------------------------------------
@@ -1870,6 +1921,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/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