Phasefield allen cahn

AUTHORS.txt

 List of contributors
 ====================
 
+Cameron Stewart
 Christian Feichtinger
 Christian Godenschwager
 Christoph Rettinger
@@ -12,14 +13,20 @@ Dominik Bartuschat
 Ehsan Fattahi
 Felix Winterhalter
 Florian Schornbaum
+Helen Schottenhamml
 Jan Götz
+Jan Hönig
+João Victor Tozatti Risso
 Johannes Habich
 Klaus Iglberger
 Kristina Pickl
 Lorenz Hufnagel
+Lukas Werner
+Markus Holzer
 Martin Bauer
 Matthias Markl
 Michael Kuron
+Nils Kohl
 Paulo Carvalho
 Regina Ammer
 Sagar Dolas
@@ -27,7 +34,9 @@ Sebastian Eibl
 Silke Bergler
 Simon Bogner
 Stefan Donath
+Stephan Seitz
 Sunil Kontham
+Tobias Leemann
 Tobias Preclik
 Tobias Scharpff
 Tobias Schruff

apps/benchmarks/CMakeLists.txt

@@ -14,9 +14,11 @@ add_subdirectory( PoiseuilleChannel )
 add_subdirectory( ProbeVsExtraMessage )
 add_subdirectory( SchaeferTurek )
 add_subdirectory( UniformGrid )
-if ( WALBERLA_BUILD_WITH_CODEGEN )
+if ( WALBERLA_BUILD_WITH_CODEGEN AND NOT WALBERLA_BUILD_WITH_CUDA )
 add_subdirectory( UniformGridGenerated )
+add_subdirectory( PhaseFieldAllenCahn )
 endif()
 if ( WALBERLA_BUILD_WITH_CUDA )
 add_subdirectory( UniformGridGPU )
-endif()
\ No newline at end of file
+add_subdirectory( PhaseFieldAllenCahn )
+endif()

apps/benchmarks/PhaseFieldAllenCahn/CMakeLists.txt

+waLBerla_link_files_to_builddir(*.prm)
+waLBerla_link_files_to_builddir(*.py)
+
+if (WALBERLA_BUILD_WITH_CUDA)
+    waLBerla_generate_target_from_python(NAME BenchmarkPhaseFieldCodeGenGPU
+            FILE multiphase_codegen.py
+            OUT_FILES initialize_phase_field_distributions.cu initialize_phase_field_distributions.h
+            initialize_velocity_based_distributions.cu initialize_velocity_based_distributions.h
+            phase_field_LB_step.cu phase_field_LB_step.h
+            hydro_LB_step.cu hydro_LB_step.h
+            PackInfo_phase_field_distributions.cu PackInfo_phase_field_distributions.h
+            PackInfo_phase_field.cu PackInfo_phase_field.h
+            PackInfo_velocity_based_distributions.cu PackInfo_velocity_based_distributions.h
+            GenDefines.h)
+
+    waLBerla_add_executable(NAME benchmark_multiphase
+            FILES benchmark_multiphase.cpp InitializerFunctions.cpp multiphase_codegen.py
+            DEPENDS blockforest core cuda field postprocessing lbm geometry timeloop gui BenchmarkPhaseFieldCodeGenGPU)
+else ()
+    waLBerla_generate_target_from_python(NAME BenchmarkPhaseFieldCodeGenCPU
+            FILE multiphase_codegen.py
+            OUT_FILES initialize_phase_field_distributions.cpp initialize_phase_field_distributions.h
+            initialize_velocity_based_distributions.cpp initialize_velocity_based_distributions.h
+            phase_field_LB_step.cpp phase_field_LB_step.h
+            hydro_LB_step.cpp hydro_LB_step.h
+            PackInfo_phase_field_distributions.cpp PackInfo_phase_field_distributions.h
+            PackInfo_phase_field.cpp PackInfo_phase_field.h
+            PackInfo_velocity_based_distributions.cpp PackInfo_velocity_based_distributions.h
+            GenDefines.h)
+
+    waLBerla_add_executable(NAME benchmark_multiphase
+            FILES benchmark_multiphase.cpp InitializerFunctions.cpp multiphase_codegen.py
+            DEPENDS blockforest core field postprocessing lbm geometry timeloop gui BenchmarkPhaseFieldCodeGenCPU)
+endif (WALBERLA_BUILD_WITH_CUDA)
+

apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file InitializerFunctions.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "python_coupling/DictWrapper.h"
+
+namespace walberla
+{
+using PhaseField_T = GhostLayerField< real_t, 1 >;
+
+void initPhaseField_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+                           const real_t R                         = 10,
+                           const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0),
+                           const real_t W                         = 5)
+{
+   for (auto& block : *blocks)
+   {
+      auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+         real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
+                          (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
+                          (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
+         phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - R) / W);
+      )
+      // clang-format on
+   }
+}
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+                        const real_t W = 5)
+{
+   auto X              = blocks->getDomainCellBB().xMax();
+   auto halfY          = (blocks->getDomainCellBB().yMax()) / 2.0;
+   double perturbation = 0.05;
+
+   for (auto& block : *blocks)
+   {
+      auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+      // clang-format off
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+         blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+         real_t tmp = perturbation * X * (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
+         phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));
+      )
+      // clang-format on
+   }
+}
+} // namespace walberla

apps/benchmarks/PhaseFieldAllenCahn/InitializerFunctions.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file InitializerFunctions.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "python_coupling/DictWrapper.h"
+#pragma once
+
+namespace walberla
+{
+void initPhaseField_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+                           Vector3< real_t > bubbleMidPoint, real_t W = 5);
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5);
+
+} // namespace walberla

apps/benchmarks/PhaseFieldAllenCahn/benchmark.py

+import os
+import waLBerla as wlb
+import pandas as pd
+
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+import sys
+
+
+class Scenario:
+    def __init__(self, timeStepStrategy, overlappingWidth, cuda_block_size):
+        # output frequencies
+        self.vtkWriteFrequency = 0
+
+        # simulation parameters
+        self.timesteps = 101
+        edge_size = 32
+        self.cells = (edge_size, edge_size, edge_size)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 1, 1)
+        self.size = (self.cells[0] * self.blocks[0],
+                     self.cells[1] * self.blocks[1],
+                     self.cells[2] * self.blocks[2])
+
+        self.timeStepStrategy = timeStepStrategy
+        self.overlappingWidth = overlappingWidth
+        self.cuda_block_size = cuda_block_size
+        self.warmupSteps = 10
+
+        # bubble parameters
+        self.bubbleRadius = min(self.size) // 4
+        self.bubbleMidPoint = (self.size[0] / 2, self.size[1] / 2, self.size[2] / 2)
+
+        self.scenario = 1  # 1 rising bubble, 2 RTI
+        self.config_dict = self.config()
+
+        self.csv_file = "benchmark.csv"
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'useGui': 0,
+                'remainingTimeLoggerFrequency': -1,
+                'timeStepStrategy': self.timeStepStrategy,
+                'overlappingWidth': self.overlappingWidth,
+                'gpuBlockSize': self.cuda_block_size,
+                'warmupSteps': self.warmupSteps,
+                'scenario': self.scenario,
+            },
+            'Boundaries_GPU': {
+                'Border': []
+            },
+            'Boundaries_CPU': {
+                'Border': []
+            },
+            'Bubble': {
+                'bubbleMidPoint': self.bubbleMidPoint,
+                'bubbleRadius': self.bubbleRadius,
+            },
+        }
+
+    @wlb.member_callback
+    def results_callback(self, **kwargs):
+        data = {}
+        data.update(self.config_dict['Parameters'])
+        data.update(self.config_dict['DomainSetup'])
+        data.update(kwargs)
+        data['executable'] = sys.argv[0]
+        data['compile_flags'] = wlb.build_info.compiler_flags
+        data['walberla_version'] = wlb.build_info.version
+        data['build_machine'] = wlb.build_info.build_machine
+        sequenceValuesToScalars(data)
+
+        df = pd.DataFrame.from_records([data])
+        if not os.path.isfile(self.csv_file):
+            df.to_csv(self.csv_file, index=False)
+        else:
+            df.to_csv(self.csv_file, index=False, mode='a', header=False)
+
+
+def overlap_benchmark():
+    scenarios = wlb.ScenarioManager()
+    overlappingWidths = [(1, 1, 1), (4, 1, 1), (8, 1, 1), (16, 1, 1), (32, 1, 1),
+                         (4, 4, 1), (8, 8, 1), (16, 16, 1), (32, 32, 1),
+                         (4, 4, 4), (8, 8, 8), (16, 16, 16), (32, 32, 32)]
+
+    cuda_blocks = [(32, 1, 1), (64, 1, 1), (128, 1, 1), (256, 1, 1), (512, 1, 1),
+                   (32, 2, 1), (64, 2, 1), (128, 2, 1), (256, 2, 1),
+                   (32, 4, 1), (64, 4, 1), (128, 4, 1),
+                   (32, 4, 2), (64, 4, 2), (128, 2, 2),
+                   (32, 8, 1), (64, 8, 1), (64, 4, 2),
+                   (32, 16, 1), (16, 16, 1), (16, 16, 2)]
+
+    # no overlap
+    scenarios.add(Scenario(timeStepStrategy='normal', overlappingWidth=(1, 1, 1), cuda_block_size=(16, 16, 1)))
+
+    # overlap
+    for overlap_strategy in ['overlap']:
+        for overlappingWidth in overlappingWidths:
+            for cuda_block in cuda_blocks:
+                scenario = Scenario(timeStepStrategy=overlap_strategy,
+                                    overlappingWidth=overlappingWidth,
+                                    cuda_block_size=cuda_block)
+                scenarios.add(scenario)
+
+
+def kernel_benchmark():
+    scenarios = wlb.ScenarioManager()
+
+    # overlap
+    # for overlap_strategy in ['phase_only', 'hydro_only', 'kernel_only']:
+    for overlap_strategy in ['overlap']:
+        scenario = Scenario(timeStepStrategy=overlap_strategy,
+                            overlappingWidth=(1, 1, 1),
+                            cuda_block_size=(128, 1, 1))
+        scenarios.add(scenario)
+
+
+# overlap_benchmark()
+kernel_benchmark()

apps/benchmarks/PhaseFieldAllenCahn/benchmark_cpu.py

+import os
+import waLBerla as wlb
+import pandas as pd
+
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+import sys
+
+
+class Scenario:
+    def __init__(self, timeStepStrategy, overlappingWidth):
+        # output frequencies
+        self.vtkWriteFrequency = -1
+
+        # simulation parameters
+        self.timesteps = 500
+        edge_size = 50
+        self.cells = (edge_size, edge_size, edge_size)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 1, 1)
+        self.size = (self.cells[0] * self.blocks[0],
+                     self.cells[1] * self.blocks[1],
+                     self.cells[2] * self.blocks[2])
+
+        self.timeStepStrategy = timeStepStrategy
+        self.overlappingWidth = overlappingWidth
+        self.cuda_block_size = (-1, -1, -1)
+        self.warmupSteps = 10
+
+        # bubble parameters
+        self.bubbleRadius = min(self.size) // 4
+        self.bubbleMidPoint = (self.size[0] / 2, self.size[1] / 2, self.size[2] / 2)
+
+        self.scenario = 1  # 1 rising bubble, 2 RTI
+        self.config_dict = self.config()
+
+        self.csv_file = "benchmark_cpu.csv"
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'useGui': 0,
+                'remainingTimeLoggerFrequency': 10.0,
+                'timeStepStrategy': self.timeStepStrategy,
+                'overlappingWidth': self.overlappingWidth,
+                'gpuBlockSize': self.cuda_block_size,
+                'warmupSteps': self.warmupSteps,
+                'scenario': self.scenario,
+            },
+            'Boundaries_GPU': {
+                'Border': []
+            },
+            'Boundaries_CPU': {
+                'Border': []
+            },
+            'Bubble': {
+                'bubbleMidPoint': self.bubbleMidPoint,
+                'bubbleRadius': self.bubbleRadius,
+            },
+        }
+
+    @wlb.member_callback
+    def results_callback(self, **kwargs):
+        data = {}
+        data.update(self.config_dict['Parameters'])
+        data.update(self.config_dict['DomainSetup'])
+        data.update(kwargs)
+        data['executable'] = sys.argv[0]
+        data['compile_flags'] = wlb.build_info.compiler_flags
+        data['walberla_version'] = wlb.build_info.version
+        data['build_machine'] = wlb.build_info.build_machine
+        sequenceValuesToScalars(data)
+
+        df = pd.DataFrame.from_records([data])
+        if not os.path.isfile(self.csv_file):
+            df.to_csv(self.csv_file, index=False)
+        else:
+            df.to_csv(self.csv_file, index=False, mode='a', header=False)
+
+
+def overlap_benchmark():
+    scenarios = wlb.ScenarioManager()
+    overlappingWidths = [(1, 1, 1), (4, 1, 1), (8, 1, 1), (16, 1, 1), (32, 1, 1),
+                         (4, 4, 1), (8, 8, 1), (16, 16, 1), (32, 32, 1),
+                         (4, 4, 4), (8, 8, 8), (16, 16, 16), (32, 32, 32),
+                         (64, 32, 32), (64, 64, 32), (64, 64, 64)]
+    # no overlap
+    scenarios.add(Scenario(timeStepStrategy='normal', overlappingWidth=(1, 1, 1)))
+
+    # overlap
+    for overlap_strategy in ['overlap']:
+        for overlappingWidth in overlappingWidths:
+            scenario = Scenario(timeStepStrategy=overlap_strategy,
+                                overlappingWidth=overlappingWidth)
+            scenarios.add(scenario)
+
+
+def kernel_benchmark():
+    scenarios = wlb.ScenarioManager()
+
+    # overlap
+    scenario = Scenario(timeStepStrategy='overlap',
+                        overlappingWidth=(8, 8, 8))
+    scenarios.add(scenario)
+
+
+# overlap_benchmark()
+kernel_benchmark()

apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file benchmark_multiphase.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformDirectScheme.h"
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/python/Exports.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "geometry/InitBoundaryHandling.h"
+
+#include "python_coupling/CreateConfig.h"
+#include "python_coupling/DictWrapper.h"
+#include "python_coupling/PythonCallback.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include <blockforest/communication/UniformBufferedScheme.h>
+
+#include "InitializerFunctions.h"
+
+//////////////////////////////
+// INCLUDE GENERATED FILES //
+////////////////////////////
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+#   include "cuda/AddGPUFieldToStorage.h"
+#   include "cuda/DeviceSelectMPI.h"
+#   include "cuda/HostFieldAllocator.h"
+#   include "cuda/NVTX.h"
+#   include "cuda/ParallelStreams.h"
+#   include "cuda/communication/GPUPackInfo.h"
+#   include "cuda/communication/MemcpyPackInfo.h"
+#   include "cuda/communication/UniformGPUScheme.h"
+
+#   include "GenDefines.h"
+#   include "PackInfo_phase_field.h"
+#   include "PackInfo_phase_field_distributions.h"
+#   include "PackInfo_velocity_based_distributions.h"
+#   include "hydro_LB_step.h"
+#   include "initialize_phase_field_distributions.h"
+#   include "initialize_velocity_based_distributions.h"
+#   include "phase_field_LB_step.h"
+
+#else
+#   include "GenDefines.h"
+#   include "PackInfo_phase_field.h"
+#   include "PackInfo_phase_field_distributions.h"
+#   include "PackInfo_velocity_based_distributions.h"
+#   include "hydro_LB_step.h"
+#   include "initialize_phase_field_distributions.h"
+#   include "initialize_velocity_based_distributions.h"
+#   include "phase_field_LB_step.h"
+#endif
+
+using namespace walberla;
+
+using PdfField_phase_T = GhostLayerField< real_t, Stencil_phase_T::Size >;
+using PdfField_hydro_T = GhostLayerField< real_t, Stencil_hydro_T::Size >;
+using VelocityField_T  = GhostLayerField< real_t, Stencil_hydro_T::Dimension >;
+using PhaseField_T     = GhostLayerField< real_t, 1 >;
+using flag_t           = walberla::uint8_t;
+using FlagField_T      = FlagField< flag_t >;
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+typedef cuda::GPUField< double > GPUField;
+#endif
+// using CommScheme_T = cuda::communication::UniformGPUScheme<stencil::D2Q9>;
+
+int main(int argc, char** argv)
+{
+   mpi::Environment env(argc, argv);
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+   cuda::selectDeviceBasedOnMpiRank();
+#endif
+
+   for (auto cfg = python_coupling::configBegin(argc, argv); cfg != python_coupling::configEnd(); ++cfg)
+   {
+      WALBERLA_MPI_WORLD_BARRIER()
+
+      auto config = *cfg;
+      logging::configureLogging(config);
+      shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGridFromConfig(config);
+
+      Vector3< uint_t > cellsPerBlock =
+         config->getBlock("DomainSetup").getParameter< Vector3< uint_t > >("cellsPerBlock");
+      // Reading parameters
+      auto parameters                    = config->getOneBlock("Parameters");
+      const std::string timeStepStrategy = parameters.getParameter< std::string >("timeStepStrategy", "normal");
+      const uint_t timesteps             = parameters.getParameter< uint_t >("timesteps", uint_c(50));
+      const real_t remainingTimeLoggerFrequency =
+         parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
+      const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
+      Vector3< int > overlappingWidth =
+         parameters.getParameter< Vector3< int > >("overlappingWidth", Vector3< int >(1, 1, 1));
+      const uint_t warmupSteps = parameters.getParameter< uint_t >("warmupSteps", uint_t(2));
+
+#if defined(WALBERLA_BUILD_WITH_CUDA)
+      // CPU fields
+      BlockDataID vel_field   = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
+      BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
+      // GPU fields
+      BlockDataID lb_phase_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
+         blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, 1);
+      BlockDataID lb_velocity_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
+         blocks, "lb velocity field on GPU", Stencil_hydro_T::Size, field::fzyx, 1);
+      BlockDataID vel_field_gpu =
+         cuda::addGPUFieldToStorage< VelocityField_T >(blocks, vel_field, "velocity field on GPU", true);
+      BlockDataID phase_field_gpu =
+         cuda::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "phase field on GPU", true);
+#else
+      BlockDataID lb_phase_field =
+         field::addToStorage< PdfField_phase_T >(blocks, "lb phase field", real_t(0), field::fzyx);