From 3a841def4d5cf75915d67cae83fc5a2d1d7fe853 Mon Sep 17 00:00:00 2001
From: Kyle Mandli <kyle.mandli@gmail.com>
Date: Tue, 15 Aug 2023 17:01:47 -0400
Subject: [PATCH] Add read function for friction.data files and a test

---
 src/python/geoclaw/data.py | 46 ++++++++++++++++++++++++++--
 tests/test_data.py         | 63 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+), 2 deletions(-)
 create mode 100644 tests/test_data.py

diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py
index 12739f31a..9e9951a9e 100755
--- a/src/python/geoclaw/data.py
+++ b/src/python/geoclaw/data.py
@@ -371,8 +371,6 @@ def write(self,data_source='setrun.py', out_file='dtopo.data'):
     def read(self, path, force=False):
         r"""Read a dtopography data file."""
 
-        print(self.dtopofiles)
-
         with open(os.path.abspath(path), 'r') as data_file:
 
             file_name = None
@@ -611,6 +609,7 @@ def __init__(self):
         # File support
         self.add_attribute('friction_files', [])
 
+
     def write(self, out_file='friction.data', data_source='setrun.py'):
 
         self.open_data_file(out_file, data_source)
@@ -647,6 +646,49 @@ def write(self, out_file='friction.data', data_source='setrun.py'):
         self.close_data_file()
 
 
+    def read(self, path="friction.data"):
+
+        super(FrictionData, self).read(path)
+
+        with open(os.path.abspath(path), 'r') as data_file:
+            # Skip header
+            for line in data_file:
+                if "#" != line[0]:
+                    break
+            value, tail = data_file.readline().split("=:")
+            self.variable_friction = bool(value)
+            value, tail = data_file.readline().split("=:")
+            self.friction_index = int(value)
+            data_file.readline()
+
+            value, tail = data_file.readline().split("=:")
+            num_friction_regions = int(value)
+            data_file.readline()
+
+            # Read in each region
+            self.friction_regions = []
+            for i in range(num_friction_regions):
+                # Read in region
+                region = []
+                for j in range(4):
+                    values, tail = data_file.readline().split("=:")
+                    region.append( [float(value) for value in values.split()] )
+                # Check
+                if len(region[3]) != len(region[2]) - 1:
+                    raise ValueError("Incorrect number of depths and coefficients.")
+
+                self.friction_regions.append(region)
+                data_file.readline()
+
+            # Read in files   
+            value, tail = data_file.readline().split("=:")
+            num_friction_files = int(value)
+            if num_friction_files > 0:
+                raise NotImplementedError("Friction files are not supported ",
+                                          "yet.")
+
+
+
 class MultilayerData(clawpack.clawutil.data.ClawData):
     r"""
     Multilayer SWE data object
diff --git a/tests/test_data.py b/tests/test_data.py
new file mode 100644
index 000000000..c9ed9bdba
--- /dev/null
+++ b/tests/test_data.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+"""Tests for reading and writing GeoClaw data files"""
+
+import tempfile
+import shutil
+import os
+
+import numpy
+
+import clawpack.clawutil.test as test
+import clawpack.geoclaw.data
+
+
+def test_read_friction_data():
+    r"""Test readinga and writing of FrictionData files"""
+
+    # Test Data
+    test_regions = []
+    test_regions.append([(-99.0, -70.0), (8.0, 32.0), 
+                       [numpy.infty, 0.0, -numpy.infty], 
+                       [0.030, 0.022]])
+    test_regions.append([(-98, 25.25), (-90, 30), 
+                       [numpy.infty, -10.0, -200.0, -numpy.infty],
+                       [0.030, 0.012, 0.022]])
+
+    # Create temp directory
+    temp_path = tempfile.mkdtemp()
+
+    try:
+        data_file = os.path.join(temp_path, "friction.data")
+        
+        # Create and write out data object
+        friction_data = clawpack.geoclaw.data.FrictionData()
+        friction_data.variable_friction = True
+        for i in range(2):
+            friction_data.friction_regions.append(test_regions[i])
+        friction_data.write(data_file)
+
+        # Read data object
+        read_friction_data = clawpack.geoclaw.data.FrictionData()
+        read_friction_data.read(data_file)
+
+        # Tests
+        for (i, region) in enumerate(test_regions):
+            for j in range(4):
+                assert numpy.allclose(region[j], 
+                                      read_friction_data.friction_regions[i][j])
+
+    except AssertionError as e:
+        # If the assertion failed then copy the contents of the directory
+        shutil.copytree(temp_path, os.path.join(os.getcwd(), 
+                                                "test_read_friction_data"))
+        raise e
+
+    finally:
+        shutil.rmtree(temp_path)
+
+
+if __name__ == "__main__":
+    test_read_friction_data()
+    print("All tests passed.")
\ No newline at end of file