Hydrothermal conductivity with depth function with fastscape

CMakeLists.txt

@@ -597,6 +597,47 @@ IF(ASPECT_WITH_LIBDAP)
   ENDIF()
 ENDIF()
 
+
+# Check for FastScape library and link it to ASPECT if requested
+SET(ASPECT_WITH_FASTSCAPE OFF CACHE BOOL "Whether the user wants to compile ASPECT with the landscape evolution code FastScape, or not.")
+MESSAGE(STATUS "Using ASPECT_WITH_FASTSCAPE = '${ASPECT_WITH_FASTSCAPE}'")
+IF(ASPECT_WITH_FASTSCAPE)
+  FIND_LIBRARY(FASTSCAPE NAMES fastscapelib_fortran PATHS $ENV{FASTSCAPE_DIR} ${FASTSCAPE_DIR} PATH_SUFFIXES lib NO_DEFAULT_PATH)
+  IF (FASTSCAPE)
+     MESSAGE(STATUS "FastScape library found at ${FASTSCAPE_DIR}")
+
+    # Get the fastscape source path so we can check the version.
+    FILE (STRINGS "${FASTSCAPE_DIR}/Makefile" _fastscape_makefile)
+    foreach(_line ${_fastscape_makefile})
+        if("${_line}" MATCHES "^CMAKE_SOURCE_DIR")
+          string(REPLACE "CMAKE_SOURCE_DIR = " "" FASTSCAPE_SOURCE_DIR ${_line})
+        endif()
+    endforeach()
+
+    # Now get the version from setup.py
+    MESSAGE(STATUS "Parsing '${FASTSCAPE_SOURCE_DIR}/setup.py' for version information")
+    FILE (STRINGS "${FASTSCAPE_SOURCE_DIR}/setup.py" _fastscape_info)
+    foreach(_line ${_fastscape_info})
+        string(STRIP ${_line} _line)
+        if("${_line}" MATCHES "^version")
+          string(REGEX REPLACE "[^0-9.]" "" FASTSCAPE_VERSION ${_line})
+        endif()
+    endforeach()
+
+    # Throw an error if the version is too old.
+    IF(${FASTSCAPE_VERSION} VERSION_LESS 2.8)
+      MESSAGE(FATAL_ERROR "The linked FastScape version is ${FASTSCAPE_VERSION}, however at least 2.8.0 is required.")
+    ENDIF()
+
+    FOREACH(_T ${TARGETS})
+      TARGET_LINK_LIBRARIES(${_T} ${FASTSCAPE})
+    ENDFOREACH()
+  ELSE()
+     MESSAGE(FATAL_ERROR "Trying to link with FastScape but libfastscapelib_fortran.so was not found in ${FASTSCAPE_DIR}")
+  ENDIF()
+ENDIF()
+
+
 #
 # NETCDF (c including parallel)
 #

contrib/fastscape/fastscape_rewrite_VTK_with_time.py

@@ -0,0 +1,156 @@
+## This is a script written by Thilo Wrona which rewrites FastScape VTK files
+## into VTS files that include the correct time of ASPECT solution files. This is
+## necessary because the VTK files do not include a time and cannot be viewed
+## simultaneously with ASPECT in paraview.
+
+import sys
+import os
+import glob
+import numpy as np
+
+import pandas as pd
+from xml.etree import ElementTree as ET
+
+#### import the simple module from the paraview, change path.
+sys.path.insert(1,'/path/to/paraview/lib/python3.8/site-packages/')
+from paraview.simple import *
+
+
+def get_times_pvd(filename):
+    # Read in the solution.pvd file
+    data = pd.read_csv(filename,header=5,delim_whitespace=True, names=[0, 1, 2, 3, 4])
+    data = data.to_numpy()
+
+     # Remove text characters and the final two rows which don't contain any time info.
+    times = np.zeros(len(data[:,1])-2)
+    for i in range(data.shape[0]-2):
+        temp_str = data[i,1]
+
+        # Column 1 contains a string of " "timestep=0" "
+        # Here we remove the first ten characters, and the final character.
+        temp_str = temp_str[10:]
+        temp_str = temp_str[:-1]
+
+        times[i] = float(temp_str)
+
+    return times 
+
+
+#%% Get file paths (absolute, not relative paths!)
+def absoluteFilePaths(directory):
+   for dirpath,_,filenames in os.walk(directory):
+       for f in filenames:
+           yield os.path.abspath(os.path.join(dirpath, f))
+
+# This will work if you run the script inside the ASPECT output folder,
+# otherwise you can put the absolute path.
+# This will make sure we only pick up the topography files.
+FileNames = glob.glob('./VTK/Topography*.vtk')
+FileNames = sorted(FileNames)
+
+DirPath = os.path.dirname(os.path.realpath(__file__))
+
+#%%
+
+# trace generated using paraview version 5.8.0
+#
+# To ensure correct image size when batch processing, please search 
+# for and uncomment the line `# renderView*.ViewSize = [*,*]`
+
+for File in FileNames:
+
+    #### disable automatic camera reset on 'Show'
+    paraview.simple._DisableFirstRenderCameraReset()
+
+    # create a new 'Legacy VTK Reader'
+    topography0000 = LegacyVTKReader(FileNames=File)
+
+    # get animation scene
+    animationScene1 = GetAnimationScene()
+
+    # get the time-keeper
+    timeKeeper1 = GetTimeKeeper()
+
+    # update animation scene based on data timesteps
+    animationScene1.UpdateAnimationUsingDataTimeSteps()
+
+    # get active view
+    renderView1 = GetActiveViewOrCreate('RenderView')
+    # uncomment following to set a specific view size
+    # renderView1.ViewSize = [882, 554]
+
+    # get layout
+    layout1 = GetLayout()
+
+    # show data in view
+    topography0000Display = Show(topography0000, renderView1, 'StructuredGridRepresentation')
+
+    # trace defaults for the display properties.
+    topography0000Display.Representation = 'Surface'
+
+    # reset view to fit data
+    renderView1.ResetCamera()
+
+    # get the material library
+    materialLibrary1 = GetMaterialLibrary()
+
+    # show color bar/color legend
+    topography0000Display.SetScalarBarVisibility(renderView1, True)
+
+    # update the view to ensure updated data information
+    renderView1.Update()
+
+    # get color transfer function/color map for 'H'
+    hLUT = GetColorTransferFunction('H')
+
+    # get opacity transfer function/opacity map for 'H'
+    hPWF = GetOpacityTransferFunction('H')
+
+    # save data
+    SaveData(File[:-4] + '.vts', proxy=topography0000, ChooseArraysToWrite=1,
+        PointDataArrays=['HHHHH','basement','catchment','drainage_area','erosion_rate','topography','total_erosion'])
+
+    #### saving camera placements for all active views
+
+    # current camera placement for renderView1
+    renderView1.CameraPosition = [225625.0, 25625.0, 878497.0779980461]
+    renderView1.CameraFocalPoint = [225625.0, 25625.0, 1135.7277145385742]
+    renderView1.CameraParallelScale = 227077.82689023562
+
+    #### uncomment the following to render all views
+    # RenderAllViews()
+    # alternatively, if you want to write images, you can use SaveScreenshot(...).
+
+    Delete(topography0000)
+    del topography0000
+
+#%% Get time steps
+times = get_times_pvd('solution.pvd')
+
+#%% Write pvd file
+
+root_attributes = {"type": "Collection", "version": "0.1", "ByteOrder": "LittleEndian"}
+
+root = ET.Element("VTKFile", attrib=root_attributes)
+root.text="\n"
+
+
+doc = ET.SubElement(root, "Collection")
+doc.text="\n"
+
+for n, File in enumerate(FileNames):
+    dataset_attributes = {"timestep": str(times[n]), "group": "", "group": "", "file": os.path.basename(File)[:-4]+'.vts'}
+
+    ET.SubElement(doc, "DataSet", attrib=dataset_attributes).text="\n"
+
+
+tree = ET.ElementTree(root)
+tree.write('./VTK/topography.pvd', xml_declaration=True, encoding='utf-8', method="xml")   
+
+
+
+
+
+
+
+

cookbooks/fastscape_eroding_box/doc/README.md

@@ -0,0 +1 @@
+This is a simple model of an eroding central block. The model is primarily used to test that the coupling installation is working properly, or for simple tests when making changes to the plugin to ensure everything is still working properly.

cookbooks/fastscape_eroding_box/fastscape_eroding_box.prm

@@ -0,0 +1,182 @@
+# This cookbook is used as a test to ensure that FastScape is properly working with ASPECT.
+# It consists of a 2.5 km high central mountain block that is eroded through the use of
+# hillslope diffusion and the Stream Power Law. Within ASPECT, the model is based off
+# convection-box-3d. However, it is set up in a way that nothing happens in the geodynamic model.
+
+# At the top, we define the number of space dimensions we would like to
+# work in:
+set Dimension                              = 3
+
+set Use years in output instead of seconds = true
+set End time                               = 2.5e5
+set Maximum time step                      = 50000
+set Output directory                       = eroding_box
+
+set Pressure normalization                 = no
+set Surface pressure                       = 0
+
+# A box that has an initial 60x60 km mountain block that is 2.5 km above the initial model surface.
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 125e3
+    set Y extent = 100e3
+    set Z extent = 25e3
+    set X repetitions = 5
+    set Y repetitions = 4
+  end
+
+  subsection Initial topography model
+    set Model name = function
+    subsection Function
+      set Function expression = if(y>20e3 && y<80e3 && x>32.5e3 && x<92.5e3, 2500, 0)
+    end
+  end
+end
+
+# We do not consider temperature in this setup
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = bottom, top
+  set List of model names = initial temperature
+end
+
+# X Set all boundaries except the top to freeslip.
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = bottom, front, back, left, right
+end
+
+# Here we setup the top boundary with fastscape.
+subsection Mesh deformation
+ set Mesh deformation boundary indicators = top : fastscape
+
+  subsection Fastscape
+    # As the highest resolution at the surface is 4 in the mesh refinement function, we set this the same.
+    set Maximum surface refinement level = 4
+
+    # Because we only have the same level across the entire surface we set this to zero.
+    # If we had multiple resolutions, e.g. 3 different cell sizes at the surface,
+    # this would be set to 2 (number of cell sizes at surface - 1).
+    set Surface refinement difference = 0
+
+    # We set FastScape as 1 level more resolved than ASPECT's initial global surface resolution.
+    set Additional fastscape refinement = 1
+
+    # Number of FastScape timesteps we want per ASPECT timestep.
+    set Number of fastscape timesteps per aspect timestep = 5
+
+    # Set the maximum allowable FastScape time step length. If the (ASPECT time step length)/(Number of FastScape steps per ASPECT step)
+    # is greater than this value, then the number of FastScape steps per ASPECT step is automatically doubled.
+    # This is continued until the FastScape timestep length is less than this value.
+    #
+    # The units of this parameter are either years or seconds, dependent on the global parameter that
+    # determines whether an input file uses one or the other.
+    set Maximum timestep length = 10000
+
+    # Seed number for initial topography noise
+    set Fastscape seed = 1000
+
+    # Because no boundaries are periodic, no flux is prescribed, and no advection occurs in FastScape,
+    # we do not need ghost nodes.
+    set Use ghost nodes = false
+    set Vertical exaggeration = 1
+
+    # Fix all boundaries in FastScape.
+    subsection Boundary conditions
+      set Front = 1
+      set Right  = 1
+      set Back    = 1
+      set Left   = 1
+    end
+
+    # We use both the Stream Power Law (SPL) and hillslope diffusion. The bedrock deposition
+    # variable allows deposition of sediment; however, because
+    # the sediment variables are not set, they are equal to
+    # bedrock values.
+    subsection Erosional parameters
+      set Drainage area exponent = 0.4
+      set Bedrock diffusivity = 1e-2
+      set Bedrock river incision rate = 1e-4
+      set Slope exponent = 1
+      set Bedrock deposition coefficient = 1
+
+      # A negative value indicates varied flow. 10 is steepest descent, and 1 is uniform distribution.
+      set Multi-direction slope exponent = -1
+    end
+  end
+
+  set Additional tangential mesh velocity boundary indicators = left, right, front, back
+
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0
+  end
+end
+
+# Dimensionless
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1
+    set Thermal expansion coefficient = 1
+    set Viscosity                     = 1
+  end
+end
+
+
+# We also have to specify that we want to use the Boussinesq
+# approximation (assuming the density in the temperature
+# equation to be constant, and incompressibility).
+subsection Formulation
+  set Formulation = Boussinesq approximation
+end
+
+
+# We set a maximum surface refinement level of 4 using the max/minimum refinement level.
+# This will make it so that the ASPECT surface refinement is the same as what we set
+# for our maximum surface refinement level for FastScape.
+subsection Mesh refinement
+  set Initial global refinement                = 4
+  set Initial adaptive refinement              = 1
+  set Time steps between mesh refinement       = 0
+  set Strategy                                 = minimum refinement function, maximum refinement function
+
+    subsection Maximum refinement function
+      set Coordinate system   = cartesian
+      set Variable names      = x,y,z
+      set Function expression = if(z>=21000, 4, 2)
+  end
+
+      subsection Minimum refinement function
+      set Coordinate system   = cartesian
+      set Variable names      = x,y,z
+      set Function expression = if(z>=21000, 4, 2)
+  end
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization
+
+  subsection Visualization
+    set Time between graphical output = 0
+    set List of output variables      = strain rate
+    set Interpolate output = true
+    set Write higher order output = true
+  end
+end