fluidity详解
fluidity详解
1.fluidity编译过程
1.1.femtools库调用方法
- 编译
fluidity/femtools目录下所有文件，打包为libfemtools.a静态库文件； - 通过
-lfemtools参数，并指定libfemtools.a静态库位置，即可调用 femtools 库内所有函数
2.fluidity主函数位置
fluidity可执行文件有4个，分别为：
- fluidity
- Burgers_Equation
- Hybridized_Helmholtz_Solver
- Shallow_Water
其中，Burgers_Equation、Hybridized_Helmholtz_Solver、Shallow_Water主程序源文件都在文件夹./main内，分别为./main/Burgers_Equation.F90，./main/Hybridized_Helmholtz_Solver.F90，./main/Shallow_Water.F90。
fluidity可执行文件源程序为最外层文件./main.cpp，main()函数则通过调用mainfl()函数来进行计算：
// Start fortran main
if(fl_command_line_options.count("simulation_name")){
mainfl();
}else{
usage(argv[0]);
exit(-1);
}
mainfl()源程序位置为./main/mainfl.F90，主要调用fluids()函数：
!######################################################
! Normal Fluidity Model
!######################################################
call tic(TICTOC_ID_SIMULATION)
ewrite(1, *) "Calling fluids from mainfl"
call fluids()
ewrite(1, *) "Exited fluids"
call toc(TICTOC_ID_SIMULATION)
call tictoc_report(2, TICTOC_ID_SIMULATION)
fluids()函数源程序位置为./main/Fluids.F90
编译fluidity可执行文件函数调用顺序为main.cpp =>./main/mainfl.F90 =>./main/Fluids.F90
3.fluidity数据结构
fluidity数据结构层次：
下层数据(quadrature_type、element_type、mesh_type)构成上层数据(element_type、mesh_type、scalar_field、vector_field、tensor_field)类型基本元素，最上层为fluidity使用的标量、矢量、矩阵等数据类型。
下面逐个介绍基本数据类型：
3.1.quadrature_type
type quadrature_type
!!< A data type which describes quadrature information. For most
!!< developers, quadrature can be treated as an opaque data type which
!!< will only be encountered when creating element_type variables to
!!< represent shape functions.
integer :: dim !! Dimension of the elements for which quadrature
!!< is required.
integer :: degree !! Degree of accuracy of quadrature.
integer :: vertices !! Number of vertices of the element.
integer :: ngi !! Number of quadrature points.
real, pointer :: weight(:)=>null() !! Quadrature weights.
real, pointer :: l(:,:)=>null() !! Locations of quadrature points.
character(len=0) :: name !! Fake name for reference counting.
!! Reference count to prevent memory leaks.
type(refcount_type), pointer :: refcount=>null()
integer :: family
end type quadrature_type
quadrature_type包含了基本单元信息，包括
- dim 维度
- degree 多项式阶数
- vertices 节点个数
- ngi 正交节点个数
- weight(😃 权重
- l(:😅 正交节点位置
- name
- refcount
- family
这些信息都是构成基本单元层面的。
!!< Given information about a quadrature, return a quad type encoding
!!< that quadrature.
function make_quadrature(vertices, dim, degree, ngi, family, stat) result (quad)
integer :: lfamily
integer :: wandzura_rule_idx, wandzura_rule_degree, max_wandzura_rule, wandzura_order
real, dimension(2, 3) :: wandzura_ref_tri
real, dimension(3, 3) :: wandzura_ref_map
real, dimension(:, :), allocatable :: tmp_coordinates
integer :: gi
integer :: gm_rule, gm_order, vertex
real, dimension(:, :), allocatable :: gm_ref_simplex
real, dimension(:, :), allocatable :: gm_ref_map
if (present(family)) then
lfamily = family
else
lfamily = FAMILY_COOLS
end if
family_if: if (lfamily == FAMILY_COOLS) then
下面根据lfamily取值对quad进行赋值，lfamily三个值分别为
- FAMILY_COOLS = 0
- FAMILY_WANDZURA = 1
- FAMILY_GM = 2
family_if: else if (lfamily == FAMILY_WANDZURA) then
! Make sure we're on triangles.
if (dim /= 2 .or. vertices /= 3) then
write (quadrature_error_message, '(a,i0,a)') ...
end if
! OK. First let's figure out which rule we want to use.
if (.not. present(degree)) then
write (quadrature_error_message, '(a,i0,a)') ...
end if
call wandzura_rule_num(max_wandzura_rule)
do wandzura_rule_idx=1,max_wandzura_rule
call wandzura_degree(wandzura_rule_idx, wandzura_rule_degree)
!! degree=idx*5
if (wandzura_rule_degree >= degree) exit
!! 当Wandzura精度超过指定精度后跳出循环
end do
if (wandzura_rule_degree < degree) then
!! 循环结束后Wandzura最大精度为30，指定精度不能超过30
write error message ..
end if
call wandzura_order_num(wandzura_rule_idx, wandzura_order)
!! 获得 wandzura_order(三角形单元中节点总个数) = ngi
call allocate(quad, vertices, wandzura_order, coords=3)
allocate(tmp_coordinates(2, wandzura_order))
quad%degree = wandzura_rule_degree
quad%dim = 2
call wandzura_rule(wandzura_rule_idx, wandzura_order, tmp_coordinates, quad%weight)
!! 获得 wandzura 节点坐标 tmp_coordinates；积分权重 quad%weight
wandzura_ref_tri(:, 1) = (/0, 0/)
wandzura_ref_tri(:, 2) = (/1, 0/)
wandzura_ref_tri(:, 3) = (/0, 1/)
call local_coords_matrix_positions(wandzura_ref_tri, wandzura_ref_map)
do gi=1,wandzura_order
quad%l(gi, 1:2) = tmp_coordinates(:, gi); quad%l(gi, 3) = 1.0
quad%l(gi, :) = matmul(wandzura_ref_map, quad%l(gi, :))
end do
在这之中有个重要的子函数调用，call allocate(quad, vertices, wandzura_order, coords=3)，目的就是为结构体quad申请内存空间。下面检查下子函数allocate的内容，
interface allocate
module procedure allocate_quad
end interface
......
subroutine allocate_quad(quad, vertices, ngi, coords, stat)
allocate(quad%weight(ngi), quad%l(ngi,coords), stat=lstat)
quad%vertices=vertices
quad%ngi=ngi
nullify(quad%refcount)
call addref(quad)
end subroutine allocate_quad
剩下最后一种定义quad方式：FAMILY_GM
family_if:elseif (lfamily == FAMILY_GM) then
......
family_if:else
......
family_if:end if
......
quad%family = lfamily
end function make_quadrature`
3.2.element_type
type element_type
!!< Type to encode shape and quadrature information for an element.
integer :: dim !! 2d or 3d?
integer :: loc !! Number of nodes.
integer :: ngi !! Number of gauss points.
integer :: degree !! Polynomial degree of element.
!! Shape functions: n is for the primitive function, dn is for partial derivatives, dn_s is for partial derivatives on surfaces.
!! n is loc x ngi, dn is loc x ngi x dim
!! dn_s is loc x ngi x face x dim
real, pointer :: n(:,:)=>null(), dn(:,:,:)=>null()
real, pointer :: n_s(:,:,:)=>null(), dn_s(:,:,:,:)=>null()
!! Polynomials defining shape functions and their derivatives.
type(polynomial), dimension(:,:), pointer :: spoly=>null(), dspoly=>null()
!! Link back to the node numbering used for this element.
type(ele_numbering_type), pointer :: numbering=>null()
!! Link back to the quadrature used for this element.
type(quadrature_type) :: quadrature
type(quadrature_type), pointer :: surface_quadrature=>null()
!! Pointer to the superconvergence data for this element.
type(superconvergence_type), pointer :: superconvergence=>null()
!! Pointer to constraints data for this element
type(constraints_type), pointer :: constraints=>null()
!! Reference count to prevent memory leaks.
type(refcount_type), pointer :: refcount=>null()
!! Dummy name to satisfy reference counting
character(len=0) :: name
end type element_type
相较而言element_type就复杂了一点，
自定义类型：ele_numbering_type，与 polynomial
type ele_numbering_type
! Type to record element numbering details.
! Differentiate tets from other elements.
integer :: faces, vertices, edges, boundaries
integer :: degree ! Degree of polynomials.
integer :: dimension ! 2D or 3D
integer :: nodes
integer :: type=ELEMENT_LAGRANGIAN
integer :: family
! Map local count coordinates to local number.
integer, dimension(:,:,:), pointer :: count2number
! Map local number to local count coordinates.
integer, dimension(:,:), pointer :: number2count
! Count coordinate which is held constant for each element boundary.
integer, dimension(:), pointer :: boundary_coord
! Value of that count coordinate on the element boundary.
integer, dimension(:), pointer :: boundary_val
end type ele_numbering_type
type polynomial
real, dimension(:), pointer :: coefs=>null()
integer :: degree=-1
end type polynomial
3.3.mesh_type
type mesh_type
!!< Mesh information for (among other things) fields.
integer, dimension(:), pointer :: ndglno
!! Flag for whether ndglno is allocated
logical :: wrapped=.true.
type(element_type) :: shape
integer :: elements
integer :: nodes
character(len=FIELD_NAME_LEN) :: name
!! path to options in the options tree
#ifdef DDEBUG
character(len=OPTION_PATH_LEN) :: option_path="/uninitialised_path/"
#else
character(len=OPTION_PATH_LEN) :: option_path
#endif
!! Degree of continuity of the field. 0 is for the conventional C0
!! discretisation. -1 for DG.
integer :: continuity=0
!! Reference count for mesh
type(refcount_type), pointer :: refcount=>null()
!! Mesh face information for those meshes (eg discontinuous) which need it.
type(mesh_faces), pointer :: faces=>null()
!! Information on subdomain_ mesh, for partially prognostic solves:
type(mesh_subdomain_mesh), pointer :: subdomain_mesh=>null()
type(adjacency_cache), pointer :: adj_lists => null()
!! array that for each node tells which column it is in
!! (column numbers usually correspond to a node number in a surface mesh)
integer, dimension(:), pointer :: columns => null()
!! if this mesh is extruded this array says which horizontal mesh element each element is below
integer, dimension(:), pointer :: element_columns => null()
!! A list of ids marking different parts of the mesh
!! so that initial conditions can be associated with it.
integer, dimension(:), pointer :: region_ids=>null()
!! Halo information for parallel simulations.
type(halo_type), dimension(:), pointer :: halos=>null()
type(halo_type), dimension(:), pointer :: element_halos=>null()
type(integer_set_vector), dimension(:), pointer :: colourings=>null()
!! A logical indicating if this mesh is periodic or not
!! (does not tell you how periodic it is... i.e. true if
!! any surface is periodic)
logical :: periodic=.false.
end type mesh_type
3.4.一维例子
test_1d.F90
function read_triangle_simple(filename, quad_degree, quad_ngi, no_faces, quad_family, mdim) result (field)
