examples/205_multivof: add sharpening_circle

examples/205_multivof/sharpening_circle/.gitignore

@@ -0,0 +1,28 @@
+velocity.pdf
+build
+add.conf
+b.dat
+*.pyc
+*.swp
+*.log
+*.status
+*.h5
+*.xmf
+*.vts
+*.pvd
+stat.dat
+*.pdf
+pid
+arg
+out
+job.id.last
+job.id
+base.conf
+mesh.conf
+np
+*.csv
+s*.vtk
+a.conf
+add.conf
+tl
+*.png

examples/205_multivof/sharpening_circle/Makefile

@@ -0,0 +1,6 @@
+m = 16 16 1
+bs = 16 16 1
+np = 1
+tl = 1440
+
+include $(shell ap.makesim)

examples/205_multivof/sharpening_circle/run_smooth

@@ -0,0 +1,8 @@
+#!/bin/sh
+set -eu
+d=smooth
+echo > add.conf
+make cleanrun
+(cd vis && rm -vf *.png && ./vis_vf.py ../s_*.vtk)
+(cd vis && mkdir -p "$d" && mv -v a_*.png "$d")
+(vis/plot_volume_error.py --output vis/$d/volume_error.pdf)

examples/205_multivof/sharpening_circle/run_stepwise

@@ -0,0 +1,8 @@
+#!/bin/sh
+set -eu
+d=stepwise
+echo "include stepwise.conf" > add.conf
+make cleanrun
+(cd vis && rm -vf *.png && ./vis_vf.py ../s_*.vtk)
+(cd vis && mkdir -p "$d" && mv -v a_*.png "$d")
+(vis/plot_volume_error.py --output vis/$d/volume_error.pdf)

examples/205_multivof/sharpening_circle/std.conf

@@ -0,0 +1,28 @@
+# 2d
+set int dim 2
+set int hypre_periodic_z 1
+
+set double extent 4
+
+# numerical
+set int sharpen 1
+set double sharpen_cfl 0.1
+set int enable_fluid 0
+set string bc_path inline
+
+# time
+set double dtmax 1
+set double tmax 10
+set double dump_field_dt 1
+set int dumpinit 1
+set string dumplist vf
+set int verbose_stages 0
+set int dumppolymarch 0
+
+# volume fraction
+set string init_vf radial_trapezoid
+set vect radial_trapezoid_c 2 2 0
+set double radial_trapezoid_rmin 0.5
+set double radial_trapezoid_rmax 1.5
+
+#include stepwise.conf

examples/205_multivof/sharpening_circle/stepwise.conf

@@ -0,0 +1,2 @@
+set double radial_trapezoid_rmin 1
+set double radial_trapezoid_rmax 1.000001

examples/205_multivof/sharpening_circle/vis/plot_volume_error.py

@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import argparse
+from argparse import Namespace
+import plottools
+import matplotlib.pyplot as plt
+import aphros
+import os
+
+myname = os.path.splitext(os.path.basename(__file__))[0]
+
+
+def load_stat(path):
+    u = np.genfromtxt(path, names=True)
+    res = Namespace(**{k: u[k] for k in u.dtype.names})
+    return res
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('dir', nargs='?', type=str, default='.')
+parser.add_argument('--output', type=str, default="volume_error.pdf")
+parser.add_argument('--ylim', nargs='*', type=float, default=[-1, 1])
+parser.add_argument('--yscpow', nargs='*', type=int, default=-15)
+args = parser.parse_args()
+
+plottools.apply_params(plt)
+
+stat = load_stat(os.path.join(args.dir, "stat.dat"))
+
+fig, ax = plt.subplots(figsize=(1.9,1.5))
+yscpow = args.yscpow
+if args.ylim:
+    ax.set_ylim(*args.ylim)
+ax.set_xlabel('step')
+ax.set_ylabel(r'volume error $[10^{{{:}}}]$'.format(yscpow))
+ax.plot(stat.step, stat.vol2_diff * 10 ** (-yscpow))
+plottools.savefig(fig, args.output)

examples/205_multivof/sharpening_circle/vis/vis_vf.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env pvbatch
+
+# state file generated using paraview version 5.9.0
+
+from paraview.simple import *
+paraview.simple._DisableFirstRenderCameraReset()
+
+import argparse
+import os
+import re
+import paratools
+
+
+parser = argparse.ArgumentParser(
+    description="Renders interface shapes from s_*.vtk files.")
+parser.add_argument('files', nargs='*', help="list of data files 's_*.vtk'")
+parser.add_argument('--force',
+                    action="store_true",
+                    help="overwrite existing files")
+args = parser.parse_args()
+
+renderView1 = CreateView('RenderView')
+renderView1.ViewSize = [1080, 1080]
+renderView1.OrientationAxesVisibility = 0
+renderView1.UseLight = 0
+renderView1.CameraPosition = [2, 2, 4]
+renderView1.CameraFocalPoint = [2, 2, 0]
+renderView1.CameraFocalDisk = 1.0
+renderView1.CameraParallelProjection = 1
+renderView1.CameraParallelScale = 1.5
+renderView1.Background = [1] * 3
+
+# https://github.com/OrdnanceSurvey/GeoDataViz-Toolkit/tree/master/Colours
+clhex_geo = [
+    "FF1F5B", "00CD6C", "009ADE", "AF58BA", "FFC61E", "F28522", "A0B1BA",
+    "A6761D", "E9002D", "FFAA00", "00B000"
+]
+
+
+def rgb(h):
+    return list(int(h[i:i + 2], 16) / 255. for i in (0, 2, 4))
+
+
+steps = paratools.ReplaceFilename(args.files, '{}', keep_dir=False)
+source_s = LegacyVTKReader(
+    FileNames=paratools.ReplaceFilename(args.files, 's_{}.vtk'))
+source_vf = XDMFReader(
+    FileNames=paratools.ReplaceFilename(args.files, 'vf_{}.xmf'))
+source_vf.CellArrayStatus = ['vf']
+sources_ft, timearrays = paratools.ApplyForceTime([source_s, source_vf])
+source_s, source_vf = sources_ft
+
+vfDisplay = Show(source_vf, renderView1)
+vfLUT = GetColorTransferFunction('vf')
+vfLUT.AutomaticRescaleRangeMode = 'Never'
+vfLUT.RGBPoints = [0.0] + [1] * 3 + [1.0] + rgb(clhex_geo[2])
+vfLUT.ColorSpace = 'RGB'
+vfLUT.NanColor = [1.0, 0.0, 0.0]
+vfLUT.Discretize = 0
+vfLUT.ScalarRangeInitialized = 1.0
+vfPWF = GetOpacityTransferFunction('vf')
+vfPWF.ScalarRangeInitialized = 1
+vfDisplay.Representation = 'Surface'
+vfDisplay.ColorArrayName = ['CELLS', 'vf']
+vfDisplay.LookupTable = vfLUT
+
+surfDisplay = Show(source_s, renderView1)
+surfDisplay.Representation = 'Wireframe'
+surfDisplay.AmbientColor = [0.0, 0.0, 0.0]
+surfDisplay.ColorArrayName = ['POINTS', '']
+surfDisplay.DiffuseColor = [0.0, 0.0, 0.0]
+surfDisplay.LineWidth = 8
+
+paratools.SetTimeStep(1, sources_ft, timearrays)
+
+paratools.SaveAnimation(steps,
+                        renderView1,
+                        sources_ft,
+                        timearrays,
+                        force=args.force)

src/client.cmake

@@ -1,7 +1,9 @@
 cmake_minimum_required(VERSION 3.3.0)
 
 # Allow using <PackageName>_ROOT variables
-cmake_policy(SET CMP0074 NEW)
+if (${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.12")
+  cmake_policy(SET CMP0074 NEW)
+endif()
 
 if (DEFINED ENV{APHROS_PREFIX})
   set(CMAKE_INSTALL_PREFIX $ENV{APHROS_PREFIX}

src/kernel/hydro.ipp

@@ -1189,7 +1189,7 @@ void Hydro<M>::InitStat(const MEB& eb) {
           "tu" + sl + "_outrate", "outlet volumetric rate of tracer " + sl, //
           [&ffv, &fcu = tracer_->GetVolumeFraction()[l], this, &meb = eb, l]() {
             Scal sum = 0;
-            mebc_fluid_.LoopBCond(meb, [&](auto cf, IdxCell c, auto bc) { //
+            mebc_fluid_.LoopBCond(meb, [&, this](auto cf, IdxCell c, auto bc) {
               if (m.IsInner(c)) {
                 if (bc.type == BCondFluidType::outlet ||
                     bc.type == BCondFluidType::outletpressure) {