Shallow-water equations on the cubed-sphere miniapp

Makefile

@@ -208,8 +208,9 @@ nekexamples  := $(OBJDIR)/nek-bps
 petscexamples.c := $(wildcard examples/petsc/*.c)
 petscexamples   := $(petscexamples.c:examples/petsc/%.c=$(OBJDIR)/petsc-%)
 # Fluid Dynamics Examples
-fluidsexamples.c := $(sort $(wildcard examples/fluids/*.c))
-fluidsexamples  := $(fluidsexamples.c:examples/fluids/%.c=$(OBJDIR)/fluids-%)
+fluidsdir        := navier-stokes shallow-water
+fluidsexamples.c := $(sort $(foreach dir,$(fluidsdir), $(wildcard examples/fluids/$(dir)/*.c)))
+fluidsexamples   := $(foreach dir,$(fluidsdir), $(fluidsexamples.c:examples/fluids/$(dir)/%.c=$(OBJDIR)/fluids-%))
 # Solid Mechanics Examples
 solidsexamples.c := $(sort $(wildcard examples/solids/*.c))
 solidsexamples   := $(solidsexamples.c:examples/solids/%.c=$(OBJDIR)/solids-%)
@@ -525,10 +526,15 @@ $(OBJDIR)/petsc-% : examples/petsc/%.c $(libceed) $(ceed.pc) | $$(@D)/.DIR
 	  PETSC_DIR="$(abspath $(PETSC_DIR))" OPT="$(OPT)" $*
 	mv examples/petsc/$* $@
 
-$(OBJDIR)/fluids-% : examples/fluids/%.c $(libceed) $(ceed.pc) | $$(@D)/.DIR
-	+$(MAKE) -C examples/fluids CEED_DIR=`pwd` \
+$(OBJDIR)/fluids-navierstokes : examples/fluids/navier-stokes/navierstokes.c $(libceed) $(ceed.pc) | $$(@D)/.DIR
+	+$(MAKE) -C examples/fluids/navier-stokes CEED_DIR=`pwd` \
 	  PETSC_DIR="$(abspath $(PETSC_DIR))" OPT="$(OPT)" $*
-	mv examples/fluids/$* $@
+	mv examples/fluids/navier-stokes/navierstokes $@
+
+$(OBJDIR)/fluids-shallowwater : examples/fluids/shallow-water/shallowwater.c $(libceed) $(ceed.pc) | $$(@D)/.DIR
+	+$(MAKE) -C examples/fluids/shallow-water CEED_DIR=`pwd` \
+	  PETSC_DIR="$(abspath $(PETSC_DIR))" OPT="$(OPT)" $*
+	mv examples/fluids/shallow-water/shallowwater $@
 
 $(OBJDIR)/solids-% : examples/solids/%.c $(libceed) $(ceed.pc) | $$(@D)/.DIR
 	+$(MAKE) -C examples/solids CEED_DIR=`pwd` \
@@ -602,6 +608,7 @@ ceedexamples : $(examples)
 nekexamples : $(nekexamples)
 mfemexamples : $(mfemexamples)
 petscexamples : $(petscexamples)
+fluidsexamples : $(fluidsexamples)
 
 # Benchmarks
 allbenchmarks = petsc-bps

doc/sphinx/source/references.bib

@@ -74,6 +74,15 @@ @article{papanastasiou1992outflow
   volume={14},
   journal={International Journal for Numerical Methods in Fluids},
   doi={10.1002/fld.1650140506}
+
+@article{rancic,
+  author = {Rančić, M. and Purser, R. J. and Mesinger, F.},
+  title = {A global shallow-water model using an expanded spherical cube: Gnomonic versus conformal coordinates},
+  journal = {Quarterly Journal of the Royal Meteorological Society},
+  volume = {122},
+  pages = {959-982},
+  doi = {10.1002/qj.49712253209},
+  year = {1996}
 }
 
 @article{straka1993numerical,
@@ -88,6 +97,25 @@ @article{straka1993numerical
   doi={10.1002/fld.1650170103}
 }
 
+@article{taylor2010ACA,
+  title={A compatible and conservative spectral element method on unstructured grids},
+  author={Mark A. Taylor and Aim{\'e} Fournier},
+  journal={J. Comput. Phys.},
+  year={2010},
+  volume={229},
+  pages={5879-5895}
+}
+
+@Article{ullrich,
+  author = {Ullrich, P. A.},
+  title = {A global finite-element shallow-water model supporting continuous and discontinuous elements},
+  journal = {Geoscientific Model Development},
+  volume = {7},
+  year = {2014},
+  pages = {3017--3035},
+  doi = {10.5194/gmd-7-3017-2014}
+}
+
 @article{williams2009roofline,
   title={Roofline: an insightful visual performance model for multicore architectures},
   author={Williams, Samuel and Waterman, Andrew and Patterson, David},

examples/fluids/.gitignore

@@ -1,5 +1,10 @@
-navierstokes
-*.vtr
-*.vts
-*.bin*
-*.vtu
+navier-stokes/navierstokes
+shallow-water/shallowwater
+navier-stokes/*.vtr
+navier-stokes/*.vts
+navier-stokes/*.bin*
+navier-stokes/*.vtu
+shallow-water/*.vtr
+shallow-water/*.vts
+shallow-water/*.bin*
+shallow-water/*.vtu

examples/fluids/README.md

@@ -1,6 +1,8 @@
+# libCEED: Fluid Dynamics Examples
+
 ## libCEED: Navier-Stokes Example
 
-This page provides a description of the Navier-Stokes example for the libCEED library, based on PETSc.
+This section provides a description of the Navier-Stokes example for the libCEED library, based on PETSc.
 
 The Navier-Stokes problem solves the compressible Navier-Stokes equations in three dimensions using an
 explicit time integration. The state variables are mass density, momentum density, and energy density.
@@ -10,6 +12,8 @@ with different problem definitions according to the application of interest.
 
 Build by using
 
+`cd navier-stokes`
+
 `make`
 
 and run with
@@ -267,3 +271,24 @@ or implicit formulation.
 
 The geometric factors and coordinate transformations required for the integration of the weak form
 are described in the file [`common.h`](common.h)
+
+## libCEED: Shallow-water Example
+
+This section provides a description of the shallow-water example for the libCEED library, based on PETSc.
+
+This solver computes the solution to the shallow-water equations on a cubed-sphere (i.e.,
+a tensor-product discrete sphere, obtained by projecting a cube inscribed in a sphere onto the surface
+of the sphere).
+
+The main shallow-water solver for libCEED is defined in [`shallowwater.c`](shallowwater.c).
+
+Build by using
+
+`cd shallow-water`
+
+`make`
+
+and run with
+
+`./shallowwater`
+

examples/fluids/index.rst

@@ -3,7 +3,7 @@
 Compressible Navier-Stokes mini-app
 ========================================
 
-This example is located in the subdirectory :file:`examples/fluids`. It solves
+This example is located in the subdirectory :file:`examples/fluids/navier-stokes`. It solves
 the time-dependent Navier-Stokes equations of compressible gas dynamics in a static
 Eulerian three-dimensional frame using unstructured high-order finite/spectral
 element spatial discretizations and explicit or implicit high-order time-stepping (available in
@@ -341,3 +341,105 @@ and non-penetration boundary conditions for :math:`\bm{u}`, and no-flux
 for mass and energy densities. This problem can be run with::
 
    ./navierstokes -problem density_current
+
+
+.. _example-petsc-shallow-water:
+
+Shallow-water Equations mini-app
+========================================
+
+This example is located in the subdirectory :file:`examples/fluids/shallow-water`. It solves
+the time-dependent shallow-water equations on a cubed-sphere. The geometry of the cubed-sphere and the relative coordinate transformations are explained in the section :ref:`example-petsc-area-sphere`. The discretization of most of the spatial differential terms needed in the shallow-water equations is already described in the section :ref:`example-petsc-bps-sphere`.
+
+The mathematical formulation (from :cite:`taylor2010ACA`) is given in what
+follows
+
+.. math::
+   :label: eq-swe
+
+   \begin{aligned}
+   \frac{\partial \bm u}{\partial t} &= - (\omega + f) \bm{\hat k} \times \bm u - \nabla \left( \frac{|\bm u|^2}{2} + g (h + h_s) \right) \\
+   \frac{\partial h}{\partial t} &= - \nabla \cdot (H_0 + h) \bm u \, , \\
+   \end{aligned}
+
+where quantities are expressed in spherical coordinates :math:`r^1 = \lambda` for longitude, :math:`r^2 = \theta` for latitude and :math:`r^3 = r` for the radius, with associated unit vectors :math:`\bm{ \hat \lambda}`, :math:`\bm {\hat \theta}`, and :math:`\bm{\hat k}`, respectively. We consider a unit sphere, hence, :math:`r = 1` and :math:`\partial / \partial r = 0`. In equations :math:numref:`eq-swe`, :math:`\bm u` represents  velocity field with longitudinal and latitudinal components :math:`\bm u = (u_{\lambda}, u_{\theta}) \equiv (u_1, u_2)`, and :math:`h` represents the height function of the fluid thickness; moreover, :math:`\omega = \nabla \times \bm u` is the vorticity, :math:`f` the Coriolis parameter such that :math:`f = 2 \bm{\Omega} \sin(\theta)` with :math:`\bm{\Omega}` Earth's rotation, :math:`g` represents the gravitational acceleration, :math:`h_s` the bottom surface (terrain) elevation, and :math:`H_0` the constant mean depth, such that the total fluid surface height is given by :math:`h_s + H_0 + h`. On the two-dimensional Riemannian manifold describing the sphere, we can use the coordinate indices :math:`x^s = {\alpha, \beta}`, so that we can express :math:`\bm u = u^{\alpha} \bm g_{\alpha} +  u^{\beta} \bm g_{\beta}`, where :math:`\bm g_{\alpha} = \partial / \partial \alpha` and :math:`\bm g_{\beta} = \partial / \partial \beta` are the natural basis vectors on the manifold :cite:`ullrich,rancic`.
+
+We solve equations :math:numref:`eq-swe` for the state variable :math:`\bm q = (q_1, q_2, q_3) = (\bm u, h) = (u_{\lambda}, u_{\theta}, h) \equiv (u_1, u_2, h)` with a semi-implicit time integration, for which 
+
+.. math::
+   :label: eq-swe-semi-implicit
+
+   \bm F(t, \bm q, \bm {\dot q}) = G(t, \bm q) \, ,
+
+where the time derivative :math:`\bm{\dot q}` is defined by
+
+.. math::
+  \bm{\dot q}(\bm q) = \sigma \bm q + \bm z
+
+in terms of :math:`\bm z` from prior state and :math:`\sigma > 0`,
+both of which depend on the specific time integration scheme. We split the time 
+integration scheme into implicit and explicit parts where the implicit terms is 
+given in what follows.
+
+.. math::
+   :label: eq-swe-implicit-part
+
+   \bm F(t, \bm q, \bm {\dot q}) := 
+   \left\{
+         \begin{array}{l}
+             F_1 (u_{\lambda}, u_{\theta}, h) = \frac{\partial u_{\lambda}}{\partial t} + g \frac{\partial }{\partial \alpha} \left( h + h_s \right)\\
+             F_2 (u_{\lambda}, u_{\theta}, h) = \frac{\partial u_{\theta}}{\partial t} + g \frac{\partial }{\partial \beta} \left( h + h_s \right)\\
+             F_3 (u_{\lambda}, u_{\theta}, h) = \frac{\partial h}{\partial t} + \frac{\partial }{\partial \alpha} \left( (H_0 + h) u_{\lambda} \right) + \frac{\partial }{\partial \beta} \left( (H_0 + h) u_{\theta} \right) .
+         \end{array}
+   \right.
+
+While for the explicit part we specify
+
+.. math::
+   :label: eq-swe-explicit-part
+
+   \bm G(t, \bm q) := 
+   \left\{
+         \begin{array}{l}
+             G_1 (u_{\lambda}, u_{\theta}, h) = - u_{\lambda} \frac{\partial u_{\lambda}}{\partial \alpha} - u_{\theta} \frac{\partial u_{\theta}}{\partial \alpha} + f u_{\theta}\\
+             G_2 (u_{\lambda}, u_{\theta}, h) =  - u_{\lambda} \frac{\partial u_{\lambda}}{\partial \beta} - u_{\theta} \frac{\partial u_{\theta}}{\partial \beta} - f u_{\lambda}\\
+             G_3 (u_{\lambda}, u_{\theta}, h) = 0 .
+         \end{array}
+   \right.
+
+The differentiation of the implicit part :math:numref:`eq-swe-implicit-part` with respect to the increment :math:`\delta \bm q = (\delta \bm u, \delta h)` gives the strong form of the incremental Jacobian :math:`\partial F_i / \partial \delta q_j`, with :math:`i,j = 1,2,3`
+
+.. math::
+   :label: eq-strong-swe-Jacobian
+
+   \frac{\partial \bm F}{\partial \delta \bm q} := 
+   \left(
+         \begin{array}{ccc}
+             \frac{\partial F_1}{\partial \delta u_{\lambda}} & \frac{\partial F_1}{\partial \delta u_{\theta}} & \frac{\partial F_1}{\partial \delta h}\\
+             \frac{\partial F_2}{\partial \delta u_{\lambda}} & \frac{\partial F_2}{\partial \delta u_{\theta}} & \frac{\partial F_2}{\partial \delta h}\\
+             \frac{\partial F_3}{\partial \delta u_{\lambda}} & \frac{\partial F_3}{\partial \delta u_{\theta}} & \frac{\partial F_3}{\partial \delta h}
+         \end{array}
+   \right)
+   =
+   \left(
+         \begin{array}{ccc}
+             0 & 0 & g \partial (\delta h) / \partial \alpha \\
+             0 & 0 & g \partial (\delta h) / \partial \beta \\
+             \frac{\partial ((H_0 + h) \delta u_{\lambda})}{\partial \alpha} & \frac{\partial ((H_0 + h) \delta u_{\theta})}{\partial \beta} & \frac{\partial (\delta h u_{\lambda})}{\partial \alpha} + \frac{\partial (\delta h u_{\theta})}{\partial \beta} ,
+         \end{array}
+   \right)
+
+whose integrands of the weak form are found by multiplying each entry of :math:`\partial F_i / \partial \delta q_j` by a test function :math:`\phi \in H^1(\Omega)`
+and integrating by parts
+
+.. math::
+   :label: eq-weak-swe-Jacobian
+
+   \left(
+         \begin{array}{ccc}
+             0 & 0 & - g \delta h \partial \phi / \partial \alpha \\
+             0 & 0 & - g \delta h \partial \phi / \partial \beta  \\
+             - \frac{\partial \phi}{\partial \alpha} ((H_0 + h) \delta u_{\lambda}) & - \frac{\partial \phi}{\partial \beta} ((H_0 + h) \delta u_{\theta}) & - \frac{\partial \phi}{\partial \alpha} (\delta h u_{\lambda}) - \frac{\partial \phi}{\partial \beta} (\delta h u_{\theta})
+         \end{array}
+   \right) .
+

examples/fluids/Makefile → examples/fluids/navier-stokes/Makefile

@@ -18,7 +18,7 @@ COMMON ?= ../../common.mk
 -include $(COMMON)
 
 PETSc.pc := $(PETSC_DIR)/$(PETSC_ARCH)/lib/pkgconfig/PETSc.pc
-CEED_DIR ?= ../..
+CEED_DIR ?= ../../..
 ceed.pc := $(CEED_DIR)/lib/pkgconfig/ceed.pc
 
 CC = $(call pkgconf, --variable=ccompiler $(PETSc.pc) $(ceed.pc))

examples/fluids/shallow-water/Makefile

@@ -0,0 +1,91 @@
+# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
+# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+# reserved. See files LICENSE and NOTICE for details.
+#
+# This file is part of CEED, a collection of benchmarks, miniapps, software
+# libraries and APIs for efficient high-order finite element and spectral
+# element discretizations for exascale applications. For more information and
+# source code availability see http://github.com/ceed.
+#
+# The CEED research is supported by the Exascale Computing Project (17-SC-20-SC)
+# a collaborative effort of two U.S. Department of Energy organizations (Office
+# of Science and the National Nuclear Security Administration) responsible for
+# the planning and preparation of a capable exascale ecosystem, including
+# software, applications, hardware, advanced system engineering and early
+# testbed platforms, in support of the nation's exascale computing imperative.
+
+PETSc.pc := $(PETSC_DIR)/$(PETSC_ARCH)/lib/pkgconfig/PETSc.pc
+CEED_DIR ?= ../../..
+ceed.pc := $(CEED_DIR)/lib/pkgconfig/ceed.pc
+
+CC = $(call pkgconf, --variable=ccompiler $(PETSc.pc) $(ceed.pc))
+CFLAGS = -std=c99 \
+  $(call pkgconf, --variable=cflags_extra $(PETSc.pc)) \
+  $(call pkgconf, --cflags-only-other $(PETSc.pc)) \
+  $(OPT)
+CPPFLAGS = $(call pkgconf, --cflags-only-I $(PETSc.pc) $(ceed.pc)) \
+  $(call pkgconf, --variable=cflags_dep $(PETSc.pc))
+LDFLAGS = $(call pkgconf, --libs-only-L --libs-only-other $(PETSc.pc) $(ceed.pc))
+LDFLAGS += $(patsubst -L%, $(call pkgconf, --variable=ldflag_rpath $(PETSc.pc))%, $(call pkgconf, --libs-only-L $(PETSc.pc) $(ceed.pc)))
+LDLIBS = $(call pkgconf, --libs-only-l $(PETSc.pc) $(ceed.pc)) -lm
+
+OBJDIR := build
+SRCDIR := src
+
+src.c := shallowwater.c $(sort $(wildcard $(SRCDIR)/*.c))
+src.o = $(src.c:%.c=$(OBJDIR)/%.o)
+
+all: shallowwater
+
+shallowwater: $(src.o) | $(PETSc.pc) $(ceed.pc)
+	$(call quiet,LINK.o) $(CEED_LDFLAGS) $^ $(LOADLIBES) $(LDLIBS) -o $@
+
+.SECONDEXPANSION: # to expand $$(@D)/.DIR
+%/.DIR :
+	@mkdir -p $(@D)
+	@touch $@
+
+# Output using the 216-color rules mode
+rule_file = $(notdir $(1))
+rule_path = $(patsubst %/,%,$(dir $(1)))
+last_path = $(notdir $(patsubst %/,%,$(dir $(1))))
+ansicolor = $(shell echo $(call last_path,$(1)) | cksum | cut -b1-2 | xargs -IS expr 2 \* S + 17)
+emacs_out = @printf "  %10s %s/%s\n" $(1) $(call rule_path,$(2)) $(call rule_file,$(2))
+color_out = @if [ -t 1 ]; then \
+				printf "  %10s \033[38;5;%d;1m%s\033[m/%s\n" \
+					$(1) $(call ansicolor,$(2)) \
+					$(call rule_path,$(2)) $(call rule_file,$(2)); else \
+				printf "  %10s %s\n" $(1) $(2); fi
+# if TERM=dumb, use it, otherwise switch to the term one
+output = $(if $(TERM:dumb=),$(call color_out,$1,$2),$(call emacs_out,$1,$2))
+
+# if V is set to non-nil, turn the verbose mode
+quiet = $(if $(V),$($(1)),$(call output,$1,$@);$($(1)))
+
+$(OBJDIR)/%.o : %.c | $$(@D)/.DIR
+	$(call quiet,CC) $(CPPFLAGS) $(CFLAGS) -c -o $@ $(abspath $<)
+
+# Rules for building the examples
+#%: %.c
+
+print: $(PETSc.pc) $(ceed.pc)
+	$(info CC      : $(CC))
+	$(info CFLAGS  : $(CFLAGS))
+	$(info CPPFLAGS: $(CPPFLAGS))
+	$(info LDFLAGS : $(LDFLAGS))
+	$(info LDLIBS  : $(LDLIBS))
+	$(info OPT     : $(OPT))
+	@true
+
+clean:
+	$(RM) -r $(OBJDIR) shallowwater *.vtu
+
+$(PETSc.pc):
+	$(if $(wildcard $@),,$(error \
+	  PETSc config not found. Please set PETSC_DIR and PETSC_ARCH))
+
+.PHONY: all print clean
+
+pkgconf = $(shell pkg-config $1 | sed -e 's/^"//g' -e 's/"$$//g')
+
+-include $(src.o:%.o=%.d)