 |
The Center for Subsurface Modeling |
|
|
The Center for Subsurface Modeling |
|
|
A joint project involving:
funded by the National Science Foundation (NSF)ParticipantsClint Dawson, Professor of Aerospace Engineering and Engineering Mechanics, The University of Texas at Austin clint@ices.utexas.edu Joannes Westerink, Professor of Civil Engineering and Geological Sciences, University of Notre Dame jjw@photius.ce.nd.edu Vadym Aizinger, Research Assistant, The University of Texas at Austin aizinger@ices.utexas.edu Jesse Feyen, Research Assistant, University of Notre Dame jfeyen@nd.edu ObjectivesAccurate mathematical modeling of coastal ocean circulation and transport has sigificant implications from an economic, environmental and public health perspective. Major inter-related issues include coastal inundation, navigation, sediment movement, pollutant transport and fisheries. Recent progress in coastal ocean modeling has emphasized two main themes. The first is the definition of relatively large computational domains which encompass much larger areas than the region of specific interest. The main concept being to place the open ocean boundaries far away in deep water non-resonant ocean basins. The second theme has been to strategically provide computational resolution using unstructured grids in order to maintain an approximately constant level of localized error throughout the domain. Thus deep water regions generally can be coarsely discretized, continental shelf and shelf break regions can be discretized with an intermediate level of grid resolution, while nearshore regions, bays, estuaries and inland floodplains require a very high level of resolution. This large domain/local grid refinement strategy has led to a number of computational difficulties. First of all, the range of flow regimes varies dramatically from the deep ocean portion of the domain to the shallow near shore and inland regions which includes inlets, rivers, and flood plains with surrounding levee systems. Not only are the depths dramatically different, but the force balances in the descriptive equations vary dramatically as well. Various algorithms perform very differently within these widely disparate flow regimes in terms of stability, accuracy and localized mass conservations properties. The second significant difficulty that is being encountered, is that the high level of grid resolution provided in localized high flow gradient and/or very shallow water depth regions actually degrades the stability properties of many algorithms that worked quite well with coarser discretizations, and work very well in regions with smoother solutions. The main focus of this project is to investigate dynamic, h-p adaptive strategies, based on mathematically-sound error estimates, for the numerical solution of shallow water models. The P.I.'s of this proposal have an extensive history in developing continuous Galerkin (CG) finite element methods for shallow water problems. We have recently investigated the use of discontinuous Galerkin (DG) methods for these problems. DG methods have a number of advantages, including flexilibity with respect to $hp$ adaptivity, and the ability to incorporate stable, higher-order approximations. We propose to further develop and analyze discontinuous Galerkin (DG) methods and exploit their strengths for solving shallow water problems which have been intractable or extremely difficult to solve by standard means. We also propose to thoroughly compare continuous and discontinuous Galerkin methods for some model problems, and to investigate novel multi-algorithmic approaches based on coupling the two methodologies for shallow water equations and related mathematical models. Our ultimate goal is to develop simulation methodology which can efficiently and accurately handle the wide range of flow regimes often observed in shallow water problems, by placing the computational power and appropriate numerical method at the right place at the right time, and to do this in an automatic way by means of suitable adaptive algorithms. Papers resulting from this project
ResultsAlgorithm development and analysisMuch of our research in this project has focused on the development of new algorithms and analyzing these algorithms for stability and convergence. This work has focused on methods for coupling continuous and discontinuous Galerkin finite elements for the 2D (depth-averaged) shallow water equations, and the extension of discontinuous Galerkin methods to 3D shallow water models. Stability and a priori error estimates for these algorithms have been derived, which demonstrate the accuracy of the methods we have developed. The Ph.D. theses of Jennifer Proft and Vadym Aizinger (available Feb 2004) at The University of Texas at Austin, have addressed these issues, in addition to several of the papers listed above. Implementation of coupled methodsThe coupled methods we have implemented can be characterized in two classes: 1) coupling across equations; i.e., using different numerical methods for different equations, and 2) coupling within equations; i.e., using domain decomposition to divide the domain into regions where we apply fully continuous or fully discontinuous approximating spaces. The first approach has been applied to the depth-integrated shallow water equations, where we apply the DG method to the primitive continuity equation, and the CG method to the momentum equation. In a recent paper (Dawson, Westerink, Pothina and Feyen 2004) we have compared a well-known CG based code, ADCIRC, with a DG based code, UTBEST, and the coupled DG/CG method. A standard test-problem, the quarter annular harbor, was chosen for this comparison. The domain is formed from a quarter annulus with an outer radius 500000 feet from the origin and an inner radius at 200000 feet. The bathymetric profile increases quadratically from 10 feet at the inner radius to 62.5 feet at the outer radius. Flow is driven by a M2 tidal forcing signal (period of 12.42 hours) which is applied uniformly at the outer radius as an elevation-specified boundary condition. The tidal forcing amplitude is ramped up using a hyperbolic tangent function from 0 to 5 feet over the first 5 days of a 30 day run. The flow propagates from the outer radius towards the inner, and since the Coriolis effect is not considered, results along all angles at a given radius should be uniform. Other details of the runs can be found in the paper. Model output is harmonically decomposed over the final 10 days of each run in order to perform the error analyses. The elevation and radial velocity time series are decomposed into tidal amplitude and phase information at each nodal location for the M2 constituent and its overtides. The numerical solutions are compared to a ``truth'' solution. In order to compute the error measures, the M2 constituent is separated into sine and cosine components for elevation and velocity. These components of the M2 elevation and velocity signals are used to obtain an error estimate and to quantify spurious modes. Model output generated by each solution method on an equilateral grid with the coarsest nodal spacing (50000 feet) is presented in the following figures. ADCIRC solution, UTBEST solution, DG/CG solution. The second approach, coupling within equations, has been applied to some one dimensional test problems in recent papers by Dawson and Proft. For example, consider the flow in a one-dimensional channel of length 100 meters, with water depth initially at a constant 1 meter. At the left boundary, we raise the depth of the water by .3 meters linearly over a period of 5 seconds. This induces a wave which flows from left to right across the domain. The steepness of the wave is influenced by the bottom friction parameter. The lower this parameter the steeper the wave. Here we consider the case of bottom friction of .06/sec. We show the DG solution obtained on a mesh of 100 equally spaced elements. Here we are plotting the elevation solution at every 50 time steps up to a final time of 25 seconds. We next show the CG solution obtained on the same mesh; notice there is substantial oscillation at the foot of the wave. The time step in both of these runs is .025 seconds. Next, we consider two coupling scenarios. First, we use the DG method in the first half of the domain and the CG method in the second half. As seen in this figure , the solution passes through the interface seamlessly, but at later time oscillations develop at the foot of the wave, similar to what is seen in the figure above. In the second scenario, we apply the CG method in the first half of the domain and the DG method in the second half. We see that the solution develops some oscillations initially at the foot of the wave. As it passes through the interface, the combined CG/DG solution shows some oscillation towards the head of the wave; these oscillations essentially vanish as the simulation proceeds, so that the final result looks very similar to the DG solution obtained above. DG for 3D shallow waterThe development of the discontinuous Galerkin method for 3D shallow water has been the subject of Vadym Aizinger's thesis work at UT Austin. Here we present one result of these simulations, that of supercritical flow in a constricted channel. Supercritical channel flows subject to a change in cross-section can lead to formation of shock and rarefaction waves. Here, we take one particular configuration, where we are given a channel of uniform depth whose boundary wall is constricted on both sides with a constriction angle of five degrees, resulting in a cross-wave pattern. Flow is induced through the left boundary with no flow through the top and bottom boundaries. The inlet Froude number is 2.5 This problem can be solved analytically; here we see the discontinuous analytical solution projected on our (smoothed) mesh. For the numerical simulations, the coarse 2D mesh consists of 3155 triangular elements with no particular orientation. We project these elements vertically to obtain a 3D grid with prismatic elements. The 3D mesh can be refined horizontally by subdividing every triangle in four using edge bisection, or/and vertically by subdividing into layers. We set vertical and horizontal eddy viscosities equal to zero and used a no-drag boundary condition at the bottom. The piecewise linear solution with 10 vertical layers shows excellent agreement with the analytical solution. Technology transferThe coupled DG/CG code developed by our group is currently being tested by the Department of Defense Naval Research Laboratory at Stennis Space Center for adoption as one of their operational models. A recent paper by C. Massey and C.A. Blain gives details on the testing of the model for problems of interest to DoD and comparisons to ADCIRC. This work resulted from collaboration and technology transfer between the DoD and our group. Last modified 14 Jan 2004 |