initial commit

geoclaw_lib/compile_geoclaw_solver.sh

+#!/bin/bash
+
+GEOCLAWSRCPATH=../src/solver/geoclaw
+
+ifort -w -g -c $GEOCLAWSRCPATH/riemannsolvers.f $GEOCLAWSRCPATH/c_bind_riemannsolvers.f90
\ No newline at end of file

src/solver/WavePropagation.hpp

+/**
+ * WavePropagation.hpp
+ *
+ ****
+ **** Abstract wave propagation solver for the Shallow Water Equations.
+ ****
+ *
+ *  Created on: Jan 03, 2012
+ *  Last Update: Jan 05, 2012
+ *
+ ****
+ *
+ *  Author: Alexander Breuer
+ *    Homepage: http://www5.in.tum.de/wiki/index.php/Dipl.-Math._Alexander_Breuer
+ *    E-Mail: breuera AT in.tum.de
+ */
+
+#ifndef WAVEPROPAGATIONSOLVER_HPP_
+#define WAVEPROPAGATIONSOLVER_HPP_
+
+namespace solver {
+  template <typename T> class WavePropagation;
+}
+
+/**
+ * Abstract wave propagation solver for the Shallow Water Equations.
+ *
+ * T should be double or float.
+ */
+template <typename T> class solver::WavePropagation {
+  protected:
+    //global variables
+    //! numerical definition of "dry".
+    const T dryTol;
+    //! gravity constant
+    const T g;
+    //! numerical definition of zero.
+    const T zeroTol;
+
+#if 0
+    //parameters for computeNetUpdates
+    T h[2];
+    T hu[2];
+    T b[2];
+    T u[2];
+
+#define hLeft  (h[0])
+#define hRight (h[1])
+
+#define huLeft  (hu[0])
+#define huRight (hu[1])
+
+#define bLeft  (b[0])
+#define bRight (b[1])
+
+#define uLeft  (u[0])
+#define uRight (u[1])
+
+#else
+    //edge-local variables
+    //! height on the left side of the edge (could change during execution).
+    T hLeft;
+    //! height on the right side of the edge (could change during execution).
+    T hRight;
+    //! momentum on the left side of the edge (could change during execution).
+    T huLeft;
+    //! momentum on the right side of the edge (could change during execution).
+    T huRight;
+    //! bathymetry on the left side of the edge (could change during execution).
+    T bLeft;
+    //! bathymetry on the right side of the edge (could change during execution).
+    T bRight;
+
+    //! velocity on the left side of the edge (computed by determineWetDryState).
+    T uLeft;
+    //! velocity on the right side of the edge (computed by determineWetDryState).
+    T uRight;
+#endif
+
+    /**
+     * The wet/dry state of the Riemann-problem.
+     */
+    enum WetDryState {
+      DryDry,               /**< Both cells are dry. */
+      WetWet,               /**< Both cells are wet. */
+      WetDryInundation,     /**< 1st cell: wet, 2nd cell: dry. 1st cell lies higher than the 2nd one. */
+      WetDryWall,           /**< 1st cell: wet, 2nd cell: dry. 1st cell lies lower than the 2nd one.
+                             *   Momentum is not large enough to overcome the difference. */
+      WetDryWallInundation, /**< 1st cell: wet, 2nd cell: dry. 1st cell lies lower than the 2nd one.
+                             *   Momentum is large enough to overcome the difference. */
+      DryWetInundation,     /**< 1st cell: dry, 2nd cell: wet. 1st cell lies lower than the 2nd one. */
+      DryWetWall,           /**< 1st cell: dry, 2nd cell: wet. 1st cell lies higher than the 2nd one.
+                             *   Momentum is not large enough to overcome the difference. */
+      DryWetWallInundation  /**< 1st cell: dry, 2nd cell: wet. 1st cell lies higher than the 2nd one.
+                             *   Momentum is large enough to overcome the difference. */
+    };
+
+    //! wet/dry state of our Riemann-problem (determined by determineWetDryState)
+    WetDryState wetDryState;
+
+
+    //! Determine the wet/dry-state and set local values if we have to.
+    virtual void determineWetDryState() = 0;
+
+    /**
+     * Constructor of a wave propagation solver.
+     *
+     * @param gravity gravity constant.
+     * @param dryTolerance numerical definition of "dry".
+     * @param zeroTolerance numerical definition of zero.
+     */
+    WavePropagation( T i_dryTolerance,
+                     T i_gravity,
+                     T i_zeroTolerance ):
+                       dryTol(i_dryTolerance),
+                       g(i_gravity),
+                       zeroTol(i_zeroTolerance) {};
+
+    /**
+     * Store parameters to member variables.
+     *
+     * @param i_hLeft height on the left side of the edge.
+     * @param i_hRight height on the right side of the edge.
+     * @param i_huLeft momentum on the left side of the edge.
+     * @param i_huRight momentum on the right side of the edge.
+     * @param i_bLeft bathymetry on the left side of the edge.
+     * @param i_bRight bathymetry on the right side of the edge.
+     */
+    inline void storeParameters( const T &i_hLeft,  const T &i_hRight,
+                                 const T &i_huLeft, const T &i_huRight,
+                                 const T &i_bLeft,  const T &i_bRight ) {
+      //store parameters to member variables
+      hLeft = i_hLeft;
+      hRight = i_hRight;
+
+      huLeft = i_huLeft;
+      huRight = i_huRight;
+
+      bLeft = i_bLeft;
+      bRight = i_bRight;
+    }
+
+    /**
+     * Store parameters to member variables.
+     *
+     * @param i_hLeft height on the left side of the edge.
+     * @param i_hRight height on the right side of the edge.
+     * @param i_huLeft momentum on the left side of the edge.
+     * @param i_huRight momentum on the right side of the edge.
+     * @param i_bLeft bathymetry on the left side of the edge.
+     * @param i_bRight bathymetry on the right side of the edge.
+     * @param i_uLeft velocity on the left side of the edge.
+     * @param i_uRight velocity on the right side of the edge.
+     */
+    inline void storeParameters( const T &i_hLeft,  const T &i_hRight,
+                                 const T &i_huLeft, const T &i_huRight,
+                                 const T &i_bLeft,  const T &i_bRight,
+                                 const T &i_uLeft,  const T &i_uRight) {
+      //store parameters to member variables
+      storeParameters( i_hLeft,  i_hRight,
+                       i_huLeft, i_huRight,
+                       i_bLeft,  i_bRight );
+
+      uLeft = i_uLeft;
+      uRight = i_uRight;
+    }
+
+  public:
+    /**
+     * Compute net updates for the cell on the left/right side of the edge.
+     * This is the default method every standalone wave propagation solver should provide.
+     *
+     * @param i_hLeft height on the left side of the edge.
+     * @param i_hRight height on the right side of the edge.
+     * @param i_huLeft momentum on the left side of the edge.
+     * @param i_huRight momentum on the right side of the edge.
+     * @param i_bLeft bathymetry on the left side of the edge.
+     * @param i_bRight bathymetry on the right side of the edge.
+     *
+     * @param o_hUpdateLeft will be set to: Net-update for the height of the cell on the left side of the edge.
+     * @param o_hUpdateRight will be set to: Net-update for the height of the cell on the right side of the edge.
+     * @param o_huUpdateLeft will be set to: Net-update for the momentum of the cell on the left side of the edge.
+     * @param o_huUpdateRight will be set to: Net-update for the momentum of the cell on the right side of the edge.
+     * @param o_maxWaveSpeed will be set to: Maximum (linearized) wave speed -> Should be used in the CFL-condition.
+     */
+    virtual void computeNetUpdates ( const T &i_hLeft,  const T &i_hRight,
+                                     const T &i_huLeft, const T &i_huRight,
+                                     const T &i_bLeft,  const T &i_bRight,
+
+                                     T &o_hUpdateLeft,
+                                     T &o_hUpdateRight,
+                                     T &o_huUpdateLeft,
+                                     T &o_huUpdateRight,
+                                     T &o_maxWaveSpeed ) = 0;
+
+    virtual ~WavePropagation() {};
+
+
+#undef hLeft
+#undef hRight
+
+#undef huLeft
+#undef huRight
+
+#undef bLeft
+#undef bRight
+
+#undef uLeft
+#undef uRight
+};
+
+#endif /* WAVEPROPAGATIONSOLVER_HPP_ */

src/unit_tests/ComponentsTest.h

+/**
+ * ComponentsTest.h
+ *
+ ****
+ **** Collection of elementary unit tests for the Augmented Riemann solver.
+ ****
+ *
+ *  Created on: Nov 30, 2011
+ *  Last Update: June 12, 2012
+ *
+ ****
+ *
+ *  Author: Alexander Breuer
+ *    Homepage: http://www5.in.tum.de/wiki/index.php/Dipl.-Math._Alexander_Breuer
+ *    E-Mail: breuera AT in.tum.de
+ *
+ ****
+ */
+
+#ifndef COMPONENTS_HPP_
+#define COMPONENTS_HPP_
+
+#define private public
+#define protected public
+#include "../solver/AugRie.hpp"
+#undef private
+#undef protected
+
+#include <string>
+
+class ComponentsTest {
+  //private:
+  //! gravity constant.
+  const double gravity;
+  //! numerical definition of "dry".
+  const double dryTol;
+  //! numerical defition of zero.
+  const double zeroTol;
+  //! maximum water height within the tests.
+  const double maxWaterHeight;
+  //! maximum wave speed within the tests.
+  //850km/h = 14166.66 m/s
+  const double maxWaveSpeed;
+  //! maximum variation of the speed within the tests.
+  //50km/h = 833.33 m/s
+  const double maxWaveSpeedVariation;
+  //! how many of the tests with random numbers should run.
+  const int numberOfRandomTests;
+  //! number of precomputed values available for the tests of the function phi.
+  const int testPhiSize;
+  //! number of phi-values, which correspond to "walls", focus: Inundation.
+  const int testPhiWallSamples;
+  //! number of phi-values, which correspond to random values.
+  const int testPhiRandomSamples;
+  //! maximum number of newton iterations in a regular setup.
+  const int maxNumberOfNewtonIterationsDefault;
+  //! maximum number of newton iterations in an expensive setup, focus: test for convergency.
+  const int maxNumberOfNewtonIterationsExpensive;
+  //! newton tolerance (when to stop iterating).
+  const double newtonTol;
+
+  //! path of the test file, which contains precomputed values.
+  const std::string testFileName;
+
+  //! Augmented Riemann solver.
+  solver::AugRie<double> mySolver;
+
+  //! water heights and velocities of the Riemann problem.
+  double hLeft, hRight, uLeft, uRight;
+
+  //! array of precomputed test values.
+  double** testPhi;
+
+  //functions
+  double createRandomNumber(const double min, const double max);
+  void readTestPhi();
+  std::string printSolverVariables();
+  std::string printSolverVariablesWithResults( const double hUpdateLeft, const double hUpdateRight,
+                                               const double huUpdateLeft, const double huUpdateRight,
+                                               const double maxWaveSpeed );
+
+
+  public:
+  ComponentsTest();
+  ~ComponentsTest();
+
+  //functions
+  void operator() ();
+
+  void checkDetermineRiemannStructure();
+  void checkComputeMiddleState();
+  void checkDetermineWetDryStateSimple();
+  void checkInundationLimitsInOneDirection(const bool i_leftCellDry);
+  void checkInundationLimits();
+  void checkComputeNetUpdatesForWallsInOneDirection(const bool i_leftCellDry);
+  void checkComputeNetUpdatesForWalls();
+  void checkComputeNetUpdatesForPositivity();
+};
+
+
+#endif /* COMPONENTS_HPP_ */

src/unit_tests/UnitTests.cpp

+/**
+ * UnitTests.cpp
+ *
+ ****
+ **** Collection of unit tests for the solvers of the shallow water equations.
+ ****
+ *
+ *  Created on: Nov 30, 2011
+ *  Last Update: June 12, 2012
+ *
+ ****
+ *
+ *  Author: Alexander Breuer
+ *    Homepage: http://www5.in.tum.de/wiki/index.php/Dipl.-Math._Alexander_Breuer
+ *    E-Mail: breuera AT in.tum.de
+ *
+ ****
+ */
+
+#include "cute.h"
+#include "ide_listener.h"
+#include "cute_runner.h"
+
+#include "ComponentsTest.h"
+
+#include "../solver/AugRieGeoClaw.hpp"
+
+#include <iostream>
+
+/**
+ * Runs the test suite of unit tests.
+ */
+void runSuite(){
+	cute::suite testSuite;
+
+	testSuite.push_back(ComponentsTest());
+	cute::ide_listener lis;
+
+	cute::makeRunner(lis)(testSuite, "running the test suite");
+}
+
+/**
+ * Main function.
+ *
+ * @return 0 if completed, 1 else.
+ */
+int main(){
+  runSuite();
+
+
+
+  /**
+   * TODO: Benchmarks will be available with the new 1D-code again (end of 2012?)
+   */
+//  benchmarks::SingleWaveOnASimpleBeach( benchmarks::OneDimensional::AugRie,
+//                                        2000,
+//                                        480,
+//                                        "/work/breuera/workspace/tsunami-benchmarks/analytical/single_wave_on_simple_beach/output/SingleWaveOnASimpleBeachAugRie.nc"
+//                                      );
+//
+//  benchmarks::SingleWaveOnASimpleBeach( benchmarks::OneDimensional::Hybrid,
+//                                        2000,
+//                                        480,
+//                                        "/work/breuera/workspace/tsunami-benchmarks/analytical/single_wave_on_simple_beach/output/SingleWaveOnASimpleBeachHybrid.nc"
+//                                       );
+//  benchmarks::DamBreakDry damBreakDry("/work/breuera/workspace/tsunami-benchmarks/analytical/single_wave_on_simple_beach/output/DamBreakDry");
+
+  std::cout << "\nUnit tests finished" << std::endl;
+
+  return 0;
+}
+
+
+

unit_tests/scripts/calculate_hstar.m

+function[hStar] = calculate_hstar(hLow, hHigh, huLow, huHigh)
+%%% constants
+zeroTol = 0.00000001;
+dryTol = 0.001;
+maxWaterHeight = 10000^3;
+
+g = 9.81;
+
+%%% Function phi
+
+phi = @(h) (h <= hLow) .* 2 .* ( sqrt(g.*h) - sqrt(g.*hLow) ) ...
+    +      (h  > hLow) .* (h - hLow) .* sqrt( g/2 .* ( 1 ./ h + 1 ./ hLow ) ) ...
+    ...
+    +      (h <= hHigh) .* 2 .* ( sqrt(g.*h) - sqrt(g.*hHigh) ) ...
+    +      (h  > hHigh) .* (h - hHigh) .* sqrt( g/2 .* ( 1 ./ h + 1 ./ hHigh ) ) ...
+    ...
+    + huHigh ./ hHigh - huLow ./ hLow;
+
+%%%
+hMin = min(hLow, hHigh);
+%hMax = max(hLow, hHigh)
+if  phi(hMin) >= 0 && ...
+    huLow / hLow - huHigh / hHigh + 2*(sqrt(g.*hLow)+sqrt(g.*hHigh)) < zeroTol
+    hStar = 0.; %large rarefaction
+    phi(hMin);
+    return;
+end
+%tic
+hStar = fzero( phi, [zeroTol, maxWaterHeight]);
+%toc
+
+%h = zeroTol: 0.01:max(hMax, hStar)+10;
+
+%plot(h, phi(h), hMin, phi(hMin), '-o', hMax, phi(hMax), '-o', hStar, 0, '-o')
+
+%text(hMin, phi(hMin), '  h_{min}')
+%text(hMax, phi(hMax), '  h_{max}')
+%text(hStar, 0, '  h_{star}')
+
+%xlabel('h')
+%ylabel('\phi(h)')
+
+%grid on
\ No newline at end of file

unit_tests/scripts/read_nc_testphi.m

+function[hLowVec, hHighVec, huLowVec, huHighVec] = read_nc_testphi(path)
+