こんばんは@soybeanです。
Fortranアドベントカレンダーの記事を拝読していて、大変勉強になります。皆様、ありがとうございます。
さて、21日目はfortranとMPIでマルチフィジクスソルバーを開発する際のクラス構造設計について議論してみたいと思います。
Who am I?
農学分野のPh.D. candidateです。土質材料ー植物間の力学的相互作用について研究しています。手法としては、構造計算や数値流体力学、計算接触力学あたりを使って現象の理解と予測を行ってきました。
加えて、植物の形態が複雑であり計算不可が大きいため、PCクラスタを構築し運用しています。
Abstract
マルチフィジクス解析のプログラムを自作する場合に、プログラムの設計は悩ましい問題です。特に、メッシュ、個々の物理場、およびソルバーを、計算リソースの不要な消費を避け、スケーラブルで、可読性を失わない形でまとめるためには多大な試行錯誤が求められます。本項では、やってみなはれ精神に基づき、すべての場の離散化を有限要素法で行う場合を想定した、マルチフィジクス解析プログラムのクラス構造を作ってみました。また、MPIへの適応性も議論します。方法として、単一の初期値境界値問題に対応したFEMDomainクラスというクラス、およびFEMDomainクラス同士の弱連成解析のターゲットとなるFEMIfaceクラスを定義します。ここで、1つのFEMDomainインスタンスは、あらゆる単一の（強連成の）初期値境界値問題について、FEMで離散化され近似解を求める際に必要十分な情報を格納できます。また、FEMIfaceインスタンスはあらゆる複数の（弱連成の）初期値境界値問題について、同様に必要十分な情報を格納できます。ここに、FEMDomainクラスは、Meshクラス（節点座標、要素コネクティビティ、要素材料番号）、MaterialPropクラス（材料パラメータ）、Boundaryクラス（境界条件、初期条件）、ShapeFunctionクラス（形状関数）、ControlParameterクラス（誤差判定に使用するパラメータ、解析ステップ数）を包含（コンポジション）しています。FEMIfaceクラスは、FEMDomainへのポインター等を包含（コンポジション）しています。実際に動かすと、各ソルバーは、強連成ソルバーであればFEMDomainインスタンスを、弱連成ソルバーであればFEMIfaceインスタンスを受け取り、連立方程式へと展開し、物理場を更新したのち、更新された情報をそれぞれ上書きすることで、全体場の更新を行い、マルチフィジクス解析を完了します。あとは、より上位のプログラムで、ソルバーを呼び出す順番を記述すれば、強連成と弱連成の入り混じった解析を行えました。スケーラビリティに関しては、ソルバーを順に呼び出すループを並列化して、FEMDomainの並列読み込みと並列更新を行うことで、MPIを用いた並列計算を実装できます。加えて、大規模計算にあたっては領域分割を行い、個々のsub領域をFEMDomainに、領域のコネクティビティ情報をFEMIfaceに格納することで、MPIを利用できます。以上により、解くべきマルチフィジクス問題をFEMDomain/FEMIfaceの両クラスのインスタンス群に分解・分類することで、計算リソースの節約、スケーラビリティおよび可読性を両立できると考えました。
1. Introduction
マルチフィジクス解析のプログラムを自作する場合に、プログラムの設計は悩ましい問題です。特に、メッシュ、個々の物理場、およびソルバーを、計算リソースの不要な消費を避け、スケーラブルで、可読性を失わない形でまとめるためには多大な試行錯誤が求められます。本項では、やってみなはれ精神に基づき、すべての場の離散化を有限要素法で行う場合を想定した、マルチフィジクス解析プログラムのクラス構造を作ってみました。また、MPIへの適応性も議論します。
2. Material and Method
方法として、単一の初期値境界値問題に対応したFEMDomainクラスというクラス、
type::FEMDomain_
type(ShapeFunction_) :: ShapeFunction
type(Mesh_) :: Mesh
type(MaterialProp_) :: MaterialProp
type(Boundary_) :: Boundary
type(ControlParameter_) :: ControlPara
real(real64) :: RealTime
character*200 :: FilePath
character*200 :: FileName
character*9 :: Dtype
character*20 :: SolverType
character*200 :: Category1
character*200 :: Category2
character*200 :: Category3
integer(int32) :: timestep
contains
procedure,public :: Init => InitializeFEMDomain
procedure,public :: Delete => DeallocateFEMDomain
procedure,public :: Export => ExportFEMDomain
procedure,public :: Import => ImportFEMDomain
procedure,public :: ImportMesh => ImportMeshFEMDomain
procedure,public :: Resize => resizeFEMDomain
procedure,public :: Merge => MergeFEMDomain
procedure,public :: AddDBoundCondition => AddDBoundCondition
procedure,public :: AddNBoundCondition => AddNBoundCondition
procedure,public :: AddTBoundCondition => AddTBoundCondition
procedure,public :: AddMaterialID => AddMaterialID
procedure,public :: SetDataType => SetDataType
procedure,public :: SetSolver => SetSolver
procedure,public :: SetName => SetName
procedure,public :: SetUp => SetUpFEMDomain
procedure,public :: InitDBC => InitDBC
procedure,public :: InitNBC => InitNBC
procedure,public :: InitTBC => InitTBC
procedure,public :: AddNBC => AddNBCFEMDomain
procedure :: MeltingSkelton => MeltingSkeltonFEMDomain
procedure,public :: SetControlPara => SetControlParaFEMDomain
procedure,public :: GmshPlotMesh => GmshPlotMesh
procedure,public :: GmshPlotContour => GmshPlotContour
procedure,public :: GmshPlotVector => GmshPlotVector
procedure,public :: GmshPlotContour2D => GmshPlotContour2D
procedure,public :: GnuplotPlotContour => GnuplotPlotContour
procedure,public :: GnuplotExportStress => GnuplotExportStress
procedure,public :: getDBCVector => getDBCVectorFEMDomain
procedure,public :: move => moveFEMDomain
procedure,public :: rotate => rotateFEMDomain
procedure,public :: meshing => meshingFEMDomain
procedure,public :: convertMeshType => convertMeshTypeFEMDomain
! for debug
procedure,public :: CheckConnectivity => CheckConnedctivityFEMDomain
end type FEMDomain_
およびFEMDomainクラス同士の弱連成解析のターゲットとなるFEMIfaceクラスを定義します。
type,extends(FEMDomain_) :: FEMDomainPointer_
type(FEMDomain_),pointer :: FEMDomainp
end type
type :: FEMIface_
type(ShapeFunction_) :: ShapeFunction1,ShapeFunction2
type(Mesh_) :: Mesh1,Mesh2 ! Mesh[12]%ElemNod is a LOCAL node pointer, not a DIRECT pointer for each domains
type(MaterialProp_) :: MaterialProp
type(ControlParameter_) :: ControlPara
type(FEMDomainPointer_),allocatable :: FEMDomains(:)
real(real64),allocatable :: NTN_NodCoord(:,:)
real(real64),allocatable :: NTS_NodCoord(:,:)
real(real64),allocatable :: STS_NodCoord(:,:)
real(real64),allocatable :: NTN_Val(:,:)
real(real64),allocatable :: NTS_Val(:,:)
real(real64),allocatable :: STS_Val(:,:)
integer(int32),allocatable :: NTN_ElemNod(:,:)
integer(int32),allocatable :: NTS_ElemNod(:,:)
integer(int32),allocatable :: STS_ElemNod(:,:)
integer(int32),allocatable :: NTN_Active(:)
integer(int32),allocatable :: NTS_Active(:)
integer(int32),allocatable :: STS_Active(:)
real(real64),allocatable :: NTN_Value(:,:)
real(real64),allocatable :: NTS_Value(:,:)
real(real64),allocatable :: STS_Value(:,:)
integer(int32),allocatable :: NTS_SegmentID(:,:)
integer(int32),allocatable :: GloNodPoint1(:,:),GloNodPoint2(:,:)
integer(int32) :: DomainID1
integer(int32) :: DomainID2
integer(int32) :: DomainID3
integer(int32) :: TimeStep
integer(int32) :: NumOfImportedDomain
character*200 :: FilePathDomain1
character*200 :: FilePathDomain2
character*200 :: FilePath
character*200 :: FileNameDomain1
character*200 :: FileNameDomain2
character*200 :: FileName
character*9 :: Dtype
character*20 :: SolverType
contains
procedure :: Init => InitializeFEMIface
procedure :: setFEMDomain => setFEMDomainFEMIface
procedure :: Delete => DeallocateFEMIface
procedure :: Import => ImportFEMIface
procedure :: GetFEMIface => GetFEMIfaceFromFEMDomains
procedure :: Export => ExportFEMIface
procedure :: GmshPlotMesh => GmshPlotMeshFEMIface
procedure :: GmshPlotNTS => GmshPlotNTSFEMIface
procedure :: GetNTNelement => GetNTNelement
procedure :: GetNTSelement => GetNTSelement
procedure :: GetGlobalNodePointer => GetGlobalNodePointerNTS
procedure :: updateTimestep => updateTimestepIface
end type
ここで、1つのFEMDomainインスタンスは、あらゆる単一の（強連成の）初期値境界値問題について、FEMで離散化され近似解を求める際に必要十分な情報を格納できます。また、FEMIfaceインスタンスはあらゆる複数の（弱連成の）初期値境界値問題について、同様に必要十分な情報を格納できます。
ここに、FEMDomainクラスは、
Meshクラス（節点座標、要素コネクティビティ、要素材料番号）、
type:: Mesh_
real(8),allocatable::NodCoord(:,:)
real(8),allocatable::NodCoordInit(:,:)
integer,allocatable::ElemNod(:,:)
integer,allocatable::FacetElemNod(:,:)
integer,allocatable::NextFacets(:,:)
integer,allocatable::SurfaceLine2D(:)
integer,allocatable::ElemMat(:)
integer,allocatable::SubMeshNodFromTo(:,:)
integer,allocatable::SubMeshElemFromTo(:,:)
integer,allocatable::SubMeshSurfFromTo(:,:)
integer :: surface
!for Interfaces
integer,allocatable::GlobalNodID(:)
character*70::ElemType
character*70 ErrorMsg
contains
procedure :: Init => InitializeMesh
procedure :: Delete => DeallocateMesh
procedure :: Copy => CopyMesh
procedure :: import => importMeshObj
procedure :: ImportElemNod => ImportElemNod
procedure :: ImportNodCoord => ImportNodCoord
procedure :: ImportElemMat => ImportElemMat
procedure :: resize => resizeMeshobj
procedure :: GetFacetElement => GetFacetElement
procedure :: GetSurface => GetSurface
procedure :: GetInterface => GetInterface
procedure :: GetInterfaceElemNod => GetInterfaceElemNod
procedure :: GetBoundingBox => GetBoundingBox
procedure :: GetFacetElemInsideBox => GetFacetElemInsideBox
procedure :: GetInterSectBox => GetInterSectBox
procedure :: GetNextFacets => GetNextFacets
procedure :: MergeMesh => MergeMesh
procedure :: ExportElemNod => ExportElemNod
procedure :: ExportNodCoord => ExportNodCoord
procedure :: ExportSurface2D => ExportSurface2D
procedure :: DisplayMesh => DisplayMesh
procedure :: ShowMesh => ShowMesh
procedure :: MeltingSkelton => MeltingSkeltonMesh
procedure :: getNumOfDomain => getNumOfDomainMesh
procedure :: SortFacet => SortFacetMesh
procedure :: Meshing => MeshingMesh
procedure :: getCircumscribedCircle => getCircumscribedCircleMesh
procedure :: getCircumscribedTriangle => getCircumscribedTriangleMesh
procedure :: removeCircumscribedTriangle => removeCircumscribedTriangleMesh
procedure :: DelauneygetNewNode => DelauneygetNewNodeMesh
procedure :: DelauneygetNewTriangle => DelauneygetNewTriangleMesh
procedure :: DelauneyremoveOverlaps => DelauneyremoveOverlapsMesh
procedure :: RemoveFailedTriangle => RemoveFailedTriangleMesh
procedure :: GetElemType => GetElemTypeMesh
procedure :: convertMeshType => ConvertMeshTypeMesh
procedure :: convertTetraToHexa => convertTetraToHexaMesh
procedure :: convertTriangleToRectangular => convertTriangleToRectangularMesh
procedure :: removeOverlappedNode =>removeOverlappedNodeMesh
end type Mesh_
MaterialPropクラス（材料パラメータ）
type :: MaterialProp_
real(real64),allocatable::MatPara(:,:)
integer(int32) :: NumOfMatPara
integer(int32) :: NumOfMaterial
character*40 :: MaterialType
character*70 ErrorMsg
contains
procedure :: Init => initializeMaterial
procedure :: Delete => DeallocateMaterialProp
end type MaterialProp_
Boundaryクラス（境界条件、初期条件）、
type::Boundary_
real(real64),allocatable::DBoundVal(:,:)
real(real64),allocatable::NBoundVal(:,:)
real(real64),allocatable::TBoundVal(:,:)
real(real64),allocatable::TBoundElemGpVal(:,:,:)
real(real64),allocatable::DBoundValInc(:,:)
real(real64),allocatable::NBoundValInc(:,:)
real(real64),allocatable::TBoundValInc(:,:)
integer(int32),allocatable::DBoundNodID(:,:)
integer(int32),allocatable::NBoundNodID(:,:)
integer(int32),allocatable::TBoundNodID(:,:)
integer(int32),allocatable::TBoundElemID(:)
integer(int32),allocatable::DBoundNum(:)
integer(int32),allocatable::NBoundNum(:)
integer(int32),allocatable::TBoundNum(:)
integer(int32),allocatable::TBoundElemNum(:)
character*70 :: ErrorMsg
contains
procedure :: Init => InitializeBoundary
procedure :: Delete => DeallocateBoundary
procedure :: CheckDataType => CheckDatatypeBoundary
procedure :: RemoveOverlap => DeleteOverlapBoundary
procedure :: ImportDBound => ImportDBound
procedure :: ImportNBound => ImportNBound
procedure :: MergeDBound => MergeDBound
procedure :: MergeNBound => MergeNBound
procedure :: removeDBC => removeDBCBoundary
procedure :: removeNBC => removeNBCBoundary
procedure :: removeTBC => removeTBCBoundary
end type Boundary_
ShapeFunctionクラス（形状関数）、
type::ShapeFunction_
real(real64),allocatable::Nmat(:)
real(real64),allocatable::dNdgzi(:,:)
real(real64),allocatable::dNdgzidgzi(:,:)
real(real64),allocatable::gzi(:)
real(real64),allocatable::GaussPoint(:,:)
real(real64),allocatable::GaussIntegWei(:)
real(real64),allocatable::Jmat(:,:),JmatInv(:,:)
real(real64),allocatable::ElemCoord(:,:)
real(real64),allocatable::ElemCoord_n(:,:)
real(real64),allocatable::du(:,:)
real(real64) :: detJ
integer(int32) :: NumOfNode
integer(int32) :: NumOfOrder
integer(int32) :: NumOfDim
integer(int32) :: NumOfGp
integer(int32) :: GpID
integer(int32) :: ierr
integer(int32) :: currentGpID
character*70::ElemType
character(len=60):: ErrorMsg
contains
procedure :: init => initShapeFunction
procedure :: update => updateShapeFunction
procedure :: SetType => SetShapeFuncType
procedure :: GetAll => GetAllShapeFunc
procedure :: Deallocate => DeallocateShapeFunction
procedure :: getType => getShapeFuncType
procedure :: GetGaussPoint => GetGaussPoint
end type ShapeFunction_
ControlParameterクラス（誤差判定に使用するパラメータ、解析ステップ数）
type::ControlParameter_
real(real64) :: Tol
integer(int32) :: SimMode
integer(int32) :: ItrTol
integer(int32) :: Timestep
contains
procedure,public :: SetControlPara => SetControlPara
end type ControlParameter_
を包含（コンポジション）しています。FEMIfaceクラスは、FEMDomainへのポインター等を包含（コンポジション）しています。
3. Results and Discussion
実際に動かすプログラムは、非常にスッキリします。
program main
use SimulatorClass
implicit none ! 美しくない一文
type(Field_),target :: world
integer,parameter :: TotalStep=40
real(8),parameter :: time=1.0d0
call world%Import()
call Simulator(world,OptionalStep=TotalStep,OptionalTime=time)
call world%Export()
end program
Simulator内に格納された各ソルバーは、強連成ソルバーであればFEMDomainインスタンスを、弱連成ソルバーであればFEMIfaceインスタンスを受け取り、連立方程式へと展開し、物理場を更新したのち、更新された情報をそれぞれ上書きすることで、全体場の更新を行い、マルチフィジクス解析を完了します。
ここに、Simulatorクラスは各種ソルバーをコンポジットしているクラスです。
type :: Simulator_
type(DiffusionEq_) ,allocatable :: DiffusionEq_Array(:)
type(FiniteDeform_),allocatable :: FiniteDeform_Array(:)
type(MultiPhysics_) ,allocatable :: MultiPhysics_Array(:)
type(ContactMechanics_),allocatable :: ContactMechanics_Array(:)
contains
procedure :: Deploy => DeploySimulator
procedure :: SetTime => SetSimulatorTime
procedure :: Run => RunSimulation
procedure :: Display => DisplaySimulation
end type
ここに、Simulatorメソッドは、ソルバー群をFEMDomain/FEMIface群にデプロイし、時間を進めます。
subroutine Simulator(world,OptionalStep,OptionalTime,SolverType)
class(Field_),target,intent(inout) :: world
type(Simulator_),target :: sim
type(MPI_) :: mpidata
integer(int32),intent(in) :: OptionalStep
real(real64),intent(in) :: OptionalTime
character(*),optional,intent(in)::SolverType
integer(int32) :: i,Step,ierr
real(real64) :: time
step=OptionalStep
time=OptionalTime
call sim%Deploy(world)
call sim%SetTime(world,time,step)
do i=1,step
call sim%Run(world,i,SolverType=SolverType)
call sim%Display(world,i)
enddo
end subroutine
ここに、FieldクラスはFEMDomain群とFEMIface群をコンポジットしているクラスです。
type :: Field_
type(FEMDomain_),allocatable::FEMDomainArray(:)
type(FEMIface_) ,allocatable::FEMIfaceArray(:)
type(FieldObjName_),allocatable::FieldList(:)
integer(int32),allocatable::Timestep(:)
real(real64),allocatable::RealTime(:)
integer(int32) :: NumberOfObject,NumberOfIface
character*200 :: FolderName
character*200 :: DomainListName
character*200 :: IfaceListName
contains
procedure :: Import => ImportField
procedure :: show => showField
procedure :: Export => ExportField
procedure :: Shift => ShiftField
procedure :: linkDomainToIface => linkDomainToIfaceField
end type
ソルバー群のFieldインスタンスへのデプロイは、ポインタでポイントすることで行います。
subroutine DeploySimulator(sim,field)
type(Field_),target,intent(inout) :: field
class(Simulator_),intent(inout) :: sim
integer(int32) :: i,j,n,DomainNum,ierr
if(.not.allocated(sim%DiffusionEq_Array) ) then
call InitializeSimulator(sim,field)
endif
n=size(field%FEMDomainArray,1)
DomainNum=n
do i=1,n
if(trim(field%FEMDomainArray(i)%SolverType)=="DiffusionEq_" )then
sim%DiffusionEq_Array(i)%FEMDomain => field%FEMDomainArray(i)
elseif(trim(field%FEMDomainArray(i)%SolverType)=="FiniteDeform_" )then
sim%FiniteDeform_Array(i)%FEMDomain => field%FEMDomainArray(i)
else
print *, "DeploySolver >> Invalid Solver Name : ",trim(field%FEMDomainArray(i)%SolverType)
stop
endif
enddo
n=size(field%FEMIfaceArray,1)
do i=1,n
if(trim(field%FEMIfaceArray(i) %SolverType)=="SyncMesh_" .or. &
trim(field%FEMIfaceArray(i)%SolverType)=="syncmesh_" )then
sim%MultiPhysics_Array(i)%FEMIFace => field%FEMIfaceArray(i)
do j=1, DomainNum
if(trim(sim%MultiPhysics_Array(i)%FEMIFace%FileNameDomain1) == &
trim(field%FEMDomainArray(j)%FileName) )then
sim%MultiPhysics_Array(i)%FEMDomain1 => field%FEMDomainArray(j)
endif
if(trim(sim%MultiPhysics_Array(i)%FEMIFace%FileNameDomain2) == &
trim(field%FEMDomainArray(j)%FileName) )then
sim%MultiPhysics_Array(i)%FEMDomain2 => field%FEMDomainArray(j)
endif
enddo
elseif(trim(field%FEMIfaceArray(i) %SolverType)=="ContactMechanics_" .or. &
trim(field%FEMIfaceArray(i)%SolverType)=="contactmechanics_" )then
sim%ContactMechanics_Array(i)%FEMIFace => field%FEMIfaceArray(i)
do j=1, DomainNum
if(trim(sim%ContactMechanics_Array(i)%FEMIFace%FileNameDomain1) == &
trim(field%FEMDomainArray(j)%FileName) )then
sim%ContactMechanics_Array(i)%FEMDomain1 => field%FEMDomainArray(j)
endif
if(trim(sim%ContactMechanics_Array(i)%FEMIFace%FileNameDomain2) == &
trim(field%FEMDomainArray(j)%FileName) )then
sim%ContactMechanics_Array(i)%FEMDomain2 => field%FEMDomainArray(j)
endif
enddo
else
print *, "DeploySolver >> Invalid Solver Name : ",trim(field%FEMIfaceArray(i) %SolverType)
stop "debug"
endif
enddo
end subroutine
スケーラビリティに関しては、ソルバーを順に呼び出すループを並列化して、FEMDomainの並列読み込みと並列更新を行うことで、MPIを用いた並列計算を実装できます。加えて、大規模計算にあたっては領域分割を行い、個々のsub領域をFEMDomainに、領域のコネクティビティ情報をFEMIfaceに格納することで、MPIを利用できます。これらの詳細な実装については今後折をみて行いたいと考えています。
また、おまけとして、MPIを呼ぶためのクラスを作りました。
type :: comment
character*200 :: comment
endtype
type:: MPI_
integer(int32) :: ierr
integer(int32) :: MyRank
integer(int32) :: PeTot
integer(int32) :: Comm1
integer(int32) :: Comm2
integer(int32) :: Comm3
integer(int32) :: Comm4
integer(int32) :: Comm5
integer(int32),allocatable::Comm(:),key(:)
integer(int32),allocatable::Stack(:,:),localstack(:)
integer(int32) :: LapTimeStep
real(real64) :: stime
real(real64) :: etime
real(real64) :: laptime(1000)
type(comment) :: comments(1000)
contains
procedure :: Start => StartMPI
procedure :: Barrier => BarrierMPI
procedure, Pass :: readMPIInt
procedure, Pass :: readMPIReal
generic :: read => readMPIInt,readMPIReal
procedure, Pass :: BcastMPIInt
procedure, Pass :: BcastMPIReal
generic :: Bcast => BcastMPIInt, BcastMPIReal
procedure, Pass :: GatherMPIInt
procedure, Pass :: GatherMPIReal
generic :: Gather => GatherMPIInt, GatherMPIReal
procedure, Pass :: ScatterMPIInt
procedure, Pass :: ScatterMPIReal
generic :: Scatter => ScatterMPIInt, ScatterMPIReal
procedure, Pass :: AllGatherMPIInt
procedure, Pass :: AllGatherMPIReal
generic :: AllGather => AllGatherMPIInt, AllGatherMPIReal
procedure, Pass :: AlltoAllMPIInt
procedure, Pass :: AlltoAllMPIReal
generic :: AlltoAll => AlltoAllMPIInt, AlltoAllMPIReal
procedure, Pass :: ReduceMPIInt
procedure, Pass :: ReduceMPIReal
generic :: Reduce => ReduceMPIInt, ReduceMPIReal
procedure, Pass :: AllReduceMPIInt
procedure, Pass :: AllReduceMPIReal
generic :: AllReduce => AllReduceMPIInt, AllReduceMPIReal
procedure :: createStack => createStackMPI
procedure :: showStack => showStackMPI
procedure :: free => freeMPI
procedure :: split => splitMPI
procedure :: copy => copyMPI
procedure :: End => EndMPI
procedure :: getLapTime => getLapTimeMPI
procedure :: showLapTime => showLapTimeMPI
procedure :: GetInfo => GetMPIInfo
end type
こちらのクラスを使うと、簡単にMPIできます。
program main
use MPIClass
implicit none
type(MPI_)::mpid
call mpid%start()
print *, "My rank is :: ", mpid%MyRank
call mpid%end()
end program
4. Conclusion
以上により、解くべきマルチフィジクス問題をFEMDomain/FEMIfaceの両クラスのインスタンス群に分解・分類することで、計算リソースの節約、スケーラビリティおよび可読性を両立できると考えました。