OPTIMISING A ONE MILLION BLOCK GRID FOR A TVZ FLOW MODEL

Abstract

The recently developed multiphase geothermal flow simulator, Waiwera, is now able to handle unstructured grids. However, to reduce errors related to the grid, certain quality factors should be optimised. In applying the finite volume method, Waiwera assumes that for each pair of neighbouring tetrahedral cells the line connecting their centroids is perfectly orthogonal to their shared interface. Any deviation from orthogonality will introduce an error in the flow simulation. Our work aims to optimize a given grid so as to minimise departures from orthogonality as much as possible. Our algorithm is written in Python. It uses the MeshIO library to read and write various formats of mesh files, and the SciPy library for optimisation algorithms. Our approach is to minimise the sum of all deviations from orthogonality within the mesh. The inputs are the node coordinates and connectivity and any surface constraints (e.g., model boundaries, faults). Then, our algorithm moves the nodes around in space in an attempt to minimise deviations from orthogonality.
Several challenges were addressed in developing the optimisation algorithm. First, we require that some nodes should be constrained to move only along boundaries (model edges, topography, faults) so that geometric features honoured by the original mesh are not altered. This reduces a node’s degrees of freedom from 3 – the number of coordinate directions - to (3 – n), where n is the number of boundaries the node lies on. As a result, the algorithm can handle plane and complex boundaries (including faults and topography), intersections between planes and between a plane and a complex surface.
Second, we aimed to minimise the running time of the algorithm so that the approach is practical for full-sized meshes (up to one million blocks). This included implementing and passing analytical Jacobians to the optimiser, and the use of multiprocessing to subdivide and parallelise optimisation of mesh subsections. The mesh is divided into m sets of disjoint (non-overlapping, non-adjacent) node clusters, where