fluidity詳解

1.fluidity編譯過程

1.1.femtools庫調(diào)用方法

編譯fluidity/femtools目錄下所有文件，打包為libfemtools.a靜態(tài)庫文件；

通過-lfemtools參數(shù)，并指定libfemtools.a靜態(tài)庫位置，即可調(diào)用 femtools 庫內(nèi)所有函數(shù)

2.fluidity主函數(shù)位置

fluidity可執(zhí)行文件有4個，分別為：

fluidity

Burgers_Equation

Hybridized_Helmholtz_Solver

Shallow_Water

其中，Burgers_Equation、Hybridized_Helmholtz_Solver、Shallow_Water主程序源文件都在文件夾./main內(nèi)，分別為./main/Burgers_Equation.F90，./main/Hybridized_Helmholtz_Solver.F90，./main/Shallow_Water.F90。

fluidity可執(zhí)行文件源程序?yàn)樽钔鈱游募?/main.cpp，main()函數(shù)則通過調(diào)用mainfl()函數(shù)來進(jìn)行計(jì)算：

// Start fortran mainif(fl_command_line_options.count("simulation_name")){mainfl(); }else{usage(argv[0]);exit(-1);}

mainfl()源程序位置為./main/mainfl.F90，主要調(diào)用fluids()函數(shù)：

!######################################################! 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()函數(shù)源程序位置為./main/Fluids.F90

編譯fluidity可執(zhí)行文件函數(shù)調(diào)用順序?yàn)閙ain.cpp =>./main/mainfl.F90 =>./main/Fluids.F90

3.fluidity數(shù)據(jù)結(jié)構(gòu)

fluidity數(shù)據(jù)結(jié)構(gòu)層次：

下層數(shù)據(jù)(quadrature_type、element_type、mesh_type)構(gòu)成上層數(shù)據(jù)(element_type、mesh_type、scalar_field、vector_field、tensor_field)類型基本元素，最上層為fluidity使用的標(biāo)量、矢量、矩陣等數(shù)據(jù)類型。
下面逐個介紹基本數(shù)據(jù)類型：

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 :: familyend type quadrature_type

quadrature_type包含了基本單元信息，包括

dim 維度

degree 多項(xiàng)式階數(shù)

vertices 節(jié)點(diǎn)個數(shù)

ngi 正交節(jié)點(diǎn)個數(shù)

weight(:) 權(quán)重

l(:,:) 正交節(jié)點(diǎn)位置

name

refcount

family

這些信息都是構(gòu)成基本單元層面的。

!!< Given information about a quadrature, return a quad type encoding!!< that quadrature.function make_quadrature(vertices, dim, degree, ngi, family, stat) result (quad)integer :: lfamilyinteger :: wandzura_rule_idx, wandzura_rule_degree, max_wandzura_rule, wandzura_orderreal, dimension(2, 3) :: wandzura_ref_trireal, dimension(3, 3) :: wandzura_ref_mapreal, dimension(:, :), allocatable :: tmp_coordinatesinteger :: giinteger :: gm_rule, gm_order, vertexreal, dimension(:, :), allocatable :: gm_ref_simplexreal, dimension(:, :), allocatable :: gm_ref_mapif (present(family)) thenlfamily = familyelselfamily = FAMILY_COOLSend iffamily_if: if (lfamily == FAMILY_COOLS) then

下面根據(jù)lfamily取值對quad進(jìn)行賦值，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) thenwrite (quadrature_error_message, '(a,i0,a)') ...end if! OK. First let's figure out which rule we want to use.if (.not. present(degree)) thenwrite (quadrature_error_message, '(a,i0,a)') ...end ifcall wandzura_rule_num(max_wandzura_rule)do wandzura_rule_idx=1,max_wandzura_rulecall wandzura_degree(wandzura_rule_idx, wandzura_rule_degree) !! degree=idx*5if (wandzura_rule_degree >= degree) exit !! 當(dāng)Wandzura精度超過指定精度后跳出循環(huán)end doif (wandzura_rule_degree < degree) then !! 循環(huán)結(jié)束后Wandzura最大精度為30，指定精度不能超過30write error message ..end ifcall wandzura_order_num(wandzura_rule_idx, wandzura_order)!! 獲得 wandzura_order(三角形單元中節(jié)點(diǎn)總個數(shù)) = ngicall allocate(quad, vertices, wandzura_order, coords=3)allocate(tmp_coordinates(2, wandzura_order))quad%degree = wandzura_rule_degreequad%dim = 2call wandzura_rule(wandzura_rule_idx, wandzura_order, tmp_coordinates, quad%weight)!! 獲得 wandzura 節(jié)點(diǎn)坐標(biāo) tmp_coordinates；積分權(quán)重 quad%weightwandzura_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_orderquad%l(gi, 1:2) = tmp_coordinates(:, gi); quad%l(gi, 3) = 1.0quad%l(gi, :) = matmul(wandzura_ref_map, quad%l(gi, :))end do

在這之中有個重要的子函數(shù)調(diào)用，call allocate(quad, vertices, wandzura_order, coords=3)，目的就是為結(jié)構(gòu)體quad申請內(nèi)存空間。下面檢查下子函數(shù)allocate的內(nèi)容，

interface allocatemodule procedure allocate_quadend interface......subroutine allocate_quad(quad, vertices, ngi, coords, stat)allocate(quad%weight(ngi), quad%l(ngi,coords), stat=lstat)quad%vertices=verticesquad%ngi=nginullify(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 = lfamilyend 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) :: quadraturetype(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 elementtype(constraints_type), pointer :: constraints=>null()!! Reference count to prevent memory leaks.type(refcount_type), pointer :: refcount=>null()!! Dummy name to satisfy reference countingcharacter(len=0) :: nameend type element_type

相較而言element_type就復(fù)雜了一點(diǎn)，

自定義類型：ele_numbering_type，與 polynomial

type ele_numbering_type! Type to record element numbering details.! Differentiate tets from other elements.integer :: faces, vertices, edges, boundariesinteger :: degree ! Degree of polynomials.integer :: dimension ! 2D or 3Dinteger :: nodesinteger :: type=ELEMENT_LAGRANGIANinteger :: 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_valend type ele_numbering_type type polynomialreal, dimension(:), pointer :: coefs=>null()integer :: degree=-1end 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 allocatedlogical :: wrapped=.true.type(element_type) :: shapeinteger :: elementsinteger :: nodescharacter(len=FIELD_NAME_LEN) :: name!! path to options in the options tree #ifdef DDEBUGcharacter(len=OPTION_PATH_LEN) :: option_path="/uninitialised_path/" #elsecharacter(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 meshtype(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 belowinteger, 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)

轉(zhuǎn)載于:https://www.cnblogs.com/li12242/p/4177175.html

總結(jié)

以上是生活随笔為你收集整理的fluidity详解的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。