fluidity详解

1.fluidity编译过程

1.1.femtools库调用方法

1. 编译fluidity/femtools目录下所有文件，打包为libfemtools.a静态库文件；
2. 通过-lfemtools参数，并指定libfemtools.a静态库位置，即可调用 femtools 库内所有函数

2.fluidity主函数位置

fluidity可执行文件有4个，分别为：

1. fluidity
2. Burgers_Equation
3. Hybridized_Helmholtz_Solver
4. 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数据结构层次：

讯享网

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包含了基本单元信息，包括

1. dim 维度
2. degree 多项式阶数
3. vertices 节点个数
4. ngi 正交节点个数
5. weight(:) 权重
6. l(:,:) 正交节点位置
7. name
8. refcount
9. 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三个值分别为

1. FAMILY_COOLS = 0
2. FAMILY_WANDZURA = 1
3. 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)