Mesh adaptation to the solution

Divergence-based The error indicator based on the area-integrated magnitude of divergence,

div;i=A_ijr v_ij A_i @u_i

@x + @v_i

@y ; (122)

is used to concentrate re…nement to areas where the velocity is highly vari-able along the ‡ow lines. This indicator is relevant not only to transcritical

‡ow with bores and standing hydraulic jumps, but also to lakes. In fact, sudden depth and width variations increase the magnitude of divergence above submerged embankments, mud banks and in straits (Figure 24, right panels).

Based on the reconstruction error of the velocity This error indicator, based on the idea of Rider and Kothe (1997), estimates the local truncation error of the velocity with

vrec;i= A_i nj

nj

X

j=1

jv_j ^v_i!jj; (123) where ^v_i!j = velocity reconstructed from cell i to cell j using the limited gradients of the ‡ow state variables, ru_i [r i;rq_i]^T:

^

v_i!j ^q_i!j

^h_i!j = q_i+rq_i (r_j r_i)

[ _i+r i (r_j r_i)] z_b;j: (124) It measures the average deviation between the velocity vector in surround-ing cells and the vector extrapolated linearly from the central cell. Hence it signals nonlinearity in the solution, where linear interpolation and cell averaging has an e¤ect on accuracy.

Each indicator is weighted by the cell areaAito …nd weaker features in coarser mesh areas of the ‡ow once the stronger features have been su¢ ciently resolved.

The three error indicators are combined into a single, nondimensional adapta-tion indicator for each leaf and parent cell i:

a_i = max

a

a;i

P_{1 sa}( _a); (125)

where the subscript a 2 [rot; div; vrec]; P_q is the q-quantile of the respective a

ands_ais the sensitivity of the adaptation to a. Standardising the error indicators using this quantile is a robust and general method that accommodates outliers and makes no assumption on how is distributed, though its calculation is costlier than that based on statistical moments . To reduce the occurrence of alternating re…nement and coarsening in subsequent adaptations, the quantile is taken on leaf cells and their parents.

4.6. MESH ADAPTATION TO THE SOLUTION 63

Figure 24. Adaptation based on the circulation (a) and divergence (b).

The modeller has control on the solution adaptation by tuning s_rot, s_div, s_vrec, acoa andnmax. A low sensitivity (sa<0.1) means that the threshold for ais raised and fewer cells are ‡agged for re…nement by that indicator. On the other hand, settings_a high (>0.25) makes it likely that the re…nement will be limited byn_max, in which case the actual value of s_a does not a¤ect the ranking of the adaptation indicator a_i any more. A lower value of the threshold parameter (acoa < 0.1) inhibits the dynamic rearrangement of re…ned areas, whereas values higher than 0.3 tend to lead to a blinking mesh, i.e., to alternating re…nement and coarsening of the same area. A particular error indicator can be switched o¤ by setting its sensitivity to zero. We found the parameter valuess_a = 0.2 and a_coa = 0.1. . . 0.2 a robust choice in general.

The adaptation procedure is as follows.

Step 1 Initiate adaptation only after the local time stepping sequence has been completed and the solution is synchronous. Calculate a;i in all cells using (121), (122) and (123). The spatial derivatives of the velocity are approx-imated using the divergence theorem (58).

Step 2 Aggregate a;P of immediate parents:

a;P =X

j

a;j; (126)

where P = parent node in the quadtree; j = leaf children of P. Step 3 Normalise a;j in cells and their parent nodes with

aa;j = ^a;j

P_{1 sa}( _a); (127)

where the quantile P_{1 sa}( _a) is based on the combined distribution of a;i

and a;P. The priority for re…nement in cells and their parent nodes is thus aj = max(adiv;j; arot;j; avrec;j): (128) Step 4 Flag each cell witha_i >1ascandidate for re…nement (CRfor short). The regularisation rules that restrict cell ratios of the initial mesh also apply to the adapted mesh. Since re…nement in a given cell may trigger the re…nement of neighbours, propagate theCR‡ag to those potentially re…ned neighbours.

Flag also the parent node of ‡agged cellsCR to prevent their coarsening.

Step 5 Flag each cell as candidate for coarsening (CC for short) if a_i < a_coa in the cell, aP < acoa in its parent node and neither the cell nor the parent has been ‡agged CR (which may occur due to regularisation). The children of the coarsened parent must be all leaves created during a previous adaptive re…nement.

In a second pass, mark thoseCC parent nodes with the C ‡ag that will be actually coarsened. For this, proceed from each CC parent and propagate theC ‡ag recursively to those neighbours that would also be coarsened due to regularisation. In this recursive process, if any of the a¤ected cells has not been ‡agged CC, cancel the C ‡ag all the way to the starting CC parent.

This guarantees that regularity will be maintained after adaptation. The total number of coarsened cells (‡aggedC) is denoted by n_coa.

Step 6 In order to control the total computational e¤ort in the simulation, identify which cells can be actually re…ned so that the maximum allowed number of cells, n_max, is not exceeded in the adapted mesh. If prior to the adaptation the mesh hasncell cells, then no more than

n_ref = max (n_max+n_coa n_cell;0) (129) number of cells may be re…ned. If this limit were exceeded by re…ning allCR

‡agged cells, restrict re…nement to those n_ref cells which have the highest priority, that is, the highest a_i. To this end, process CR cells in the order of decreasing a_i and mark them with the R ‡ag by propagating the R ‡ag also to those cells that will be re…ned due to mesh regularisation, until the

4.6. MESH ADAPTATION TO THE SOLUTION 65 limit n_ref is reached or all CR cells have been ‡agged. For convenience, the limit to the total cell count (nmax) is given as a multiple of the cell count in the initial mesh.

Step 7 Calculate and store rui for cells ‡agged R for re…nement. Then split these cells into four equal-size children.

Step 8 Merge the four quadrant cells of parents ‡aggedC. Beforehand, integrate the ‡ow state variables in wettable children i and store the result for the parent P;

(Au)_P =X

i

A_iu_i: (130)

At this point the mesh has been reshaped.

The properties of new cells (bottom elevation, attributes) are not derived from the existing mesh but re-rendered from the digital elevation model and the bound-ary polygons. This allows the model to identify, as the mesh is being re…ned, small-scale bed and coverage features that were previously smeared due to insu¢ cient resolution.

The solution, on the other hand, is projected to the new mesh by linear in-terpolation based on the last u …eld. Since the bottom surface has also changed, surface elevations must be adjusted so that volume is conserved during adapta-tion. This adjustment interferes with the water surface in fully wet cells. The fact that zb is assigned by area-weighted averaging of the digital elevation model minimises the interference, however, if any quadrant was dry before coarsening or a new quadrant of an originally wet cell became dry after re…nement, then the adjustment of will perturb the solution locally. Likewise, the interpolation of p and q to the new mesh also introduces minor numerical shocks (even in fully wet cells) because momentum conservation is not enforced.

The remaining steps of the adaptation procedure that update the …eld variables are given next.

Step 9 Re-render the re…nement attributes, ‡ow attributes and bed topography to the new mesh.

Step 10 Cells that were previously ‡agged R have become parent nodes to the four new children created during re…nement. Disaggregate the u values to these children by linear extrapolation as

u_i =u_P +ru_P(r_i r_P); (131) where the gradient ru_P has been calculated and saved in Step 7.

If volume conservation is enforced during adaptation, distribute the volume h_PA_P in the children iteratively so thatP

h_iA_i =h_PA_P taking into account that some of the four children may have become dry or exterior cells, and

…nish by calculating the surface elevation, _i =hi+zb;i.

Step 11 Aggregate the solution in the coarsened cells (indexedinow butP before adaptation) using the integral variables saved in Step 8 as

u_i = (Au)_i=A_i: (132) Again, if volume conservation is required then the depth and not the surface elevation is integrated over the cell area in Step 8,

(Ah)_P =X

i

A_ih_i; (133)

and that volume is used to restore the depth of the coarsened cell:

i = (Ah)_i=A_i z_b;i: (134) Step 12 Adaptation may have altered the mesh in shallow regions where the reinterpreted bed elevation in some previously dry cells are now below the surrounding water surface. Also, near the ‡ow domain boundary some new cells whose parent used to be exterior may have obtained wettable status and the solution must be initialised there. Though setting h = 0 in these newly wettable cells conserves water volume, the small shock waves that would inevitably arise due to the discontinuity of the water surface would lead to oscillatory water boundaries. Hence it is desirable to relax exact volume conservation and extrapolate the solution to the newly ‡ooded cells.

Dry cells that have just been re…ned are recursively ‡ooded with inverse distance weighting from the surrounding wet cells as

i = X

j

j=d_ij X

j

(1=d_ij)

; q_i = (0;0)^T; (135)

where i denotes the dry cell in question; j 2 wet neighbours of i; d_ij = jr_i r_jj is the distance of the cell centres.

In document Solution-adaptive 2D modelling of wind-induced lake circulation (Pldal 79-84)