diff --git a/apps/showcases/Thermocapillary/InitializerFunctions.cpp b/apps/showcases/Thermocapillary/InitializerFunctions.cpp
index 47b147497d73ee8cdfb327068fe0c67fc2df6c14..464e56d1939a37f2524c61455e8799e8271985cd 100644
--- a/apps/showcases/Thermocapillary/InitializerFunctions.cpp
+++ b/apps/showcases/Thermocapillary/InitializerFunctions.cpp
@@ -47,15 +47,17 @@ void initPhaseFieldDroplet(const shared_ptr< StructuredBlockStorage >& blocks, B
void initMicroChannel(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, BlockDataID temperatureFieldID,
const real_t Th, const real_t T0, const real_t Tc, const real_t W)
{
- auto halfX = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
- auto halfY = real_c((blocks->getDomainCellBB().yMax())) / real_c(2.0);
+ auto domain = blocks->getDomain();
+ auto halfX = real_c(domain.xMax()) / real_c(2.0);
+ auto halfY = real_c(domain.yMax()) / real_c(2.0);
for (auto& block : *blocks)
{
auto phaseField = block.getData< ScalarField_T >(phaseFieldID);
// clang-format off
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh( (real_c(globalCell[1]) - halfY) / (W / real_c(2.0)) ));
+ auto globalCellCenter = blocks->getCellCenter(globalCell);
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh( (real_c(globalCellCenter[1]) - halfY) / (W / real_c(2.0)) ));
)
// clang-format on
@@ -65,8 +67,8 @@ void initMicroChannel(const shared_ptr< StructuredBlockStorage >& blocks, BlockD
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
if(globalCell[1] == -1)
{
- auto X = (globalCell[0] < 0) ? 0.0 : globalCell[0];
- temperatureField->get(x, y, z) = Th + T0 * cos( (math::pi / halfX) * (X - halfX) );
+ auto globalCellCenter = blocks->getCellCenter(globalCell);
+ temperatureField->get(x, y, z) = Th + T0 * cos( (math::pi / halfX) * (globalCellCenter[0] - halfX) );
}
else
{
diff --git a/apps/showcases/Thermocapillary/microchannel2D.py b/apps/showcases/Thermocapillary/microchannel2D.py
index 9b183ace47dee2573461b24eb79f177dd3874769..a710eb7d4cf80dcf98f6c2192c77f47fcb99c65c 100755
--- a/apps/showcases/Thermocapillary/microchannel2D.py
+++ b/apps/showcases/Thermocapillary/microchannel2D.py
@@ -14,11 +14,11 @@ from waLBerla.tools.sqlitedb import sequenceValuesToScalars
class Scenario:
def __init__(self, heat_solver_rk_or_lbm, case):
- self.reference_length = 256
+ self.reference_length = 80
# simulation parameters
- self.timesteps = 1
- self.pre_thermal_timesteps = 1
+ self.timesteps = 400000
+ self.pre_thermal_timesteps = 100000
self.vtkWriteFrequency = 0
self.dbWriteFrequency = 10000
@@ -28,7 +28,7 @@ class Scenario:
self.periodic = (1, 0, 0)
self.heat_solver_rk_or_lbm = heat_solver_rk_or_lbm
- self.order_rk_solver = 4
+ self.order_rk_solver = 2
self.case = case
@@ -81,7 +81,8 @@ class Scenario:
self.u_y = u_y
self.t_a = t_a
l0 = self.reference_length
- self.XX, self.YY = np.meshgrid(np.arange(self.domain_size[0]) - l0, np.arange(self.domain_size[1]) - l0 // 2)
+ self.XX, self.YY = np.meshgrid(np.arange(self.domain_size[0]) - (l0-0.5),
+ np.arange(self.domain_size[1]) - (l0 // 2 - 0.5))
self.config_dict = self.config()
@@ -197,8 +198,8 @@ class Scenario:
contour_2 = plt.contour(self.XX, self.YY, temperature_field.T, levels=levels, colors=['k'])
clabels = ax.clabel(contour_2, inline=1, fontsize=10, fmt='%2.0lf')
[txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=0)) for txt in clabels]
- plt.ylim((-128, 128))
- plt.xlim((-256, 256))
+ plt.ylim((-self.reference_length//2, self.reference_length//2))
+ plt.xlim((-self.reference_length, self.reference_length))
plt.savefig(f'{path}/temperature_{t}.png', dpi=300)
plt.close('all')
@@ -214,8 +215,8 @@ class Scenario:
velocity_field[::n, ::n, 0].T, velocity_field[::n, ::n, 1].T,
angles='xy', scale_units='xy', scale=0.00001, color='b')
plt.grid()
- plt.ylim((-128, 128))
- plt.xlim((-256, 256))
+ plt.ylim((-self.reference_length, self.reference_length))
+ plt.xlim((-2*self.reference_length, 2*self.reference_length))
plt.savefig(f'{path}/velocity_{t}.png', dpi=300)
plt.close('all')
@@ -249,4 +250,4 @@ scenarios = wlb.ScenarioManager()
# scenarios.add(Scenario(heat_solver_rk_or_lbm=False, case=1))
scenarios.add(Scenario(heat_solver_rk_or_lbm=True, case=1))
# scenarios.add(Scenario(heat_solver_rk_or_lbm=False, case=2))
-scenarios.add(Scenario(heat_solver_rk_or_lbm=True, case=2))
+# scenarios.add(Scenario(heat_solver_rk_or_lbm=True, case=2))
diff --git a/apps/showcases/Thermocapillary/microchannel3D.py b/apps/showcases/Thermocapillary/microchannel3D.py
index 0726d7114ed6cb60c759ba090e694dabff1edba8..e6b53f4b814b1d6fb46bc23692f37405ff265d56 100755
--- a/apps/showcases/Thermocapillary/microchannel3D.py
+++ b/apps/showcases/Thermocapillary/microchannel3D.py
@@ -14,7 +14,7 @@ from waLBerla.tools.sqlitedb import sequenceValuesToScalars
class Scenario:
def __init__(self, heat_solver_rk_or_lbm, case):
- self.reference_length = 256
+ self.reference_length = 80
# simulation parameters
self.timesteps = 400001
@@ -22,7 +22,7 @@ class Scenario:
self.vtkWriteFrequency = 0
self.dbWriteFrequency = 10000
- self.cells = (self.reference_length, self.reference_length, 10)
+ self.cells = (self.reference_length, self.reference_length, 3)
self.blocks = (2, 1, 1)
self.domain_size = tuple([c * b for c, b in zip(self.cells, self.blocks)])
self.periodic = (1, 0, 1)
@@ -46,7 +46,7 @@ class Scenario:
self.temperature_h = 20
self.temperature_l = 4
- self.interface_thickness = 5
+ self.interface_thickness = 4
self.mobility = 0.05
self.velocity_wall = 0
else:
@@ -64,7 +64,7 @@ class Scenario:
self.temperature_h = 20
self.temperature_l = 4
- self.interface_thickness = 5
+ self.interface_thickness = 4
self.mobility = 0.05
self.velocity_wall = 0
@@ -80,7 +80,8 @@ class Scenario:
self.u_y = u_y
self.t_a = t_a
l0 = self.reference_length
- self.XX, self.YY = np.meshgrid(np.arange(self.domain_size[0]) - l0, np.arange(self.domain_size[1]) - l0 // 2)
+ self.XX, self.YY = np.meshgrid(np.arange(self.domain_size[0]) - (l0-0.5),
+ np.arange(self.domain_size[1]) - (l0 // 2 - 0.5))
self.config_dict = self.config()
@@ -195,8 +196,8 @@ class Scenario:
contour_2 = plt.contour(self.XX, self.YY, temperature_field.T, levels=levels, colors=['k'])
clabels = ax.clabel(contour_2, inline=1, fontsize=10, fmt='%2.0lf')
[txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=0)) for txt in clabels]
- plt.ylim((-128, 128))
- plt.xlim((-256, 256))
+ plt.ylim((-self.reference_length//2, self.reference_length//2))
+ plt.xlim((-self.reference_length, self.reference_length))
plt.savefig(f'{path}/temperature_{t}.png', dpi=300)
plt.close('all')
@@ -214,8 +215,8 @@ class Scenario:
angles='xy', scale_units='xy', scale=0.00001)
ax1.set_title("lattice Boltzmann")
plt.grid()
- plt.ylim((-128, 128))
- plt.xlim((-256, 256))
+ plt.ylim((-self.reference_length, self.reference_length))
+ plt.xlim((-2*self.reference_length, 2*self.reference_length))
plt.savefig(f'{path}/velocity_{t}.png', dpi=300)
plt.close('all')
diff --git a/apps/showcases/Thermocapillary/thermocapillary.cpp b/apps/showcases/Thermocapillary/thermocapillary.cpp
index fd8b5d155b9043e0d0df492d3820ab1a3a0aa4ea..3a6611f98e3b133609682d97cf6a66e22b887881 100644
--- a/apps/showcases/Thermocapillary/thermocapillary.cpp
+++ b/apps/showcases/Thermocapillary/thermocapillary.cpp
@@ -68,11 +68,13 @@ using FlagField_T = FlagField< flag_t >;
auto DiffusionCallback = [](const Cell& pos, const shared_ptr< StructuredBlockForest >& blocks, IBlock& block,
const real_t Th, const real_t T0) {
- auto x_half = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
- Cell global_cell;
- blocks->transformBlockLocalToGlobalCell(global_cell, block, pos);
+ auto domain = blocks->getDomain();
+ auto halfX = real_c(domain.xMax()) / real_c(2.0);
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, pos);
+ auto globalCellCenter = blocks->getCellCenter(globalCell);
- return Th + T0 * cos( (math::pi / x_half) * (real_c(global_cell[0]) - x_half) );
+ return Th + T0 * cos( (math::pi / halfX) * (real_c(globalCellCenter[0]) - halfX) );
};
int main(int argc, char** argv)
@@ -549,9 +551,9 @@ int main(int argc, char** argv)
DiffusionInitialisation = std::bind(DiffusionCallback, _1, _2, _3, t_h, t_0);
#if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
- lbm::thermal_LB_DiffusionDirichlet_static thermal_LB_DiffusionStatic(blocks, temperature_PDFs_ID, vel_field_gpu, temperature_ref);
+ lbm::thermal_LB_DiffusionDirichlet_static thermal_LB_DiffusionStatic(blocks, temperature_PDFs_ID, temperature_ref);
thermal_LB_DiffusionStatic.fillFromFlagField< FlagField_T >(blocks, flag_field_thermal_LB_ID, diffusionStaticFlagUID, fluidFlagUID);
- lbm::thermal_LB_DiffusionDirichlet_dynamic thermal_LB_DiffusionDynamic(blocks, temperature_PDFs_ID, vel_field_gpu, DiffusionInitialisation);
+ lbm::thermal_LB_DiffusionDirichlet_dynamic thermal_LB_DiffusionDynamic(blocks, temperature_PDFs_ID, DiffusionInitialisation);
thermal_LB_DiffusionDynamic.fillFromFlagField< FlagField_T >(blocks, flag_field_thermal_LB_ID, diffusionDynamicFlagUID, fluidFlagUID);
lbm::thermal_LB_NeumannByCopy thermal_LB_Neumann(blocks, temperature_PDFs_ID);
thermal_LB_Neumann.fillFromFlagField< FlagField_T >(blocks, flag_field_thermal_LB_ID, neumannFlagUID, fluidFlagUID);
@@ -559,10 +561,9 @@ int main(int argc, char** argv)
pystencils::ContactAngle contact_angle(blocks, phase_field_gpu, interface_thickness, angle);
contact_angle.fillFromFlagField< FlagField_T >(blocks, flag_field_allen_cahn_LB_ID, wallFlagUID, fluidFlagUID);
#else
- lbm::thermal_LB_DiffusionDirichlet_static thermal_LB_DiffusionStatic(blocks, temperature_PDFs_ID,
- velocity_field_ID, temperature_ref);
+ lbm::thermal_LB_DiffusionDirichlet_static thermal_LB_DiffusionStatic(blocks, temperature_PDFs_ID, temperature_ref);
thermal_LB_DiffusionStatic.fillFromFlagField< FlagField_T >(blocks, flag_field_thermal_LB_ID, diffusionStaticFlagUID, fluidFlagUID);
- lbm::thermal_LB_DiffusionDirichlet_dynamic thermal_LB_DiffusionDynamic(blocks, temperature_PDFs_ID, velocity_field_ID, DiffusionInitialisation);
+ lbm::thermal_LB_DiffusionDirichlet_dynamic thermal_LB_DiffusionDynamic(blocks, temperature_PDFs_ID, DiffusionInitialisation);
thermal_LB_DiffusionDynamic.fillFromFlagField< FlagField_T >(blocks, flag_field_thermal_LB_ID, diffusionDynamicFlagUID, fluidFlagUID);
lbm::thermal_LB_NeumannByCopy thermal_LB_Neumann(blocks, temperature_PDFs_ID);
thermal_LB_Neumann.fillFromFlagField< FlagField_T >(blocks, flag_field_thermal_LB_ID, neumannFlagUID, fluidFlagUID);
diff --git a/apps/showcases/Thermocapillary/thermocapillary_codegen.py b/apps/showcases/Thermocapillary/thermocapillary_codegen.py
index 1930f0019a19c139c8432a34794553e389799a4f..58c3c9f37a51a3cfbbfc52f91c3f67aef226acef 100644
--- a/apps/showcases/Thermocapillary/thermocapillary_codegen.py
+++ b/apps/showcases/Thermocapillary/thermocapillary_codegen.py
@@ -11,7 +11,7 @@ from pystencils.simp import (add_subexpressions_for_field_reads, sympy_cse, inse
from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil, create_lb_update_rule
from lbmpy.creationfunctions import create_lb_method
-from lbmpy.boundaries import DiffusionDirichlet, NeumannByCopy, NoSlip, UBB
+from lbmpy.boundaries import FixedDensity, NeumannByCopy, NoSlip, UBB
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
@@ -23,7 +23,6 @@ from lbmpy.phasefield_allen_cahn.parameter_calculation import ThermocapillaryPar
import pystencils_walberla
from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field, generate_info_header
-from lbmpy_walberla.additional_data_handler import DiffusionDirichletAdditionalDataHandler
from lbmpy_walberla import generate_boundary, generate_lb_pack_info
# to disable type conversion warnings
@@ -44,9 +43,9 @@ with CodeGeneration() as ctx:
contact_angle_in_degrees = sp.Symbol("alpha")
generate_with_heat_source = False
- stencil_phase = LBStencil(Stencil.D3Q19)
- stencil_hydro = LBStencil(Stencil.D3Q19)
- stencil_thermal = LBStencil(Stencil.D3Q7)
+ stencil_phase = LBStencil(Stencil.D2Q9)
+ stencil_hydro = LBStencil(Stencil.D2Q9)
+ stencil_thermal = LBStencil(Stencil.D2Q9)
assert (stencil_phase.D == stencil_hydro.D == stencil_thermal.D)
full_stencil = LBStencil(Stencil.D2Q9) if stencil_phase.D == 2 else LBStencil(Stencil.D3Q27)
@@ -374,14 +373,12 @@ const bool GeneratedHeatSource = false;
varying_parameters=vp, gpu_indexing_params=sweep_params,
cpu_openmp=openmp, cpu_vectorize_info=cpu_vec)
- dirichlet_bc_static = DiffusionDirichlet(TypedSymbol('T_c', dtype=field_type), u, data_type=field_type)
+ dirichlet_bc_static = FixedDensity(TypedSymbol('T_c', dtype=field_type), u, data_type=field_type)
generate_boundary(ctx, 'thermal_LB_DiffusionDirichlet_static', dirichlet_bc_static, method_thermal,
target=target, streaming_pattern='pull', cpu_openmp=openmp, data_type=field_type_pdfs)
- dirichlet_bc_dynamic = DiffusionDirichlet(lambda *args: None, u, data_type=field_type)
- diffusion_data_handler = DiffusionDirichletAdditionalDataHandler(stencil_thermal, dirichlet_bc_dynamic)
+ dirichlet_bc_dynamic = FixedDensity(lambda *args: None, u, data_type=field_type)
generate_boundary(ctx, 'thermal_LB_DiffusionDirichlet_dynamic', dirichlet_bc_dynamic, method_thermal,
- additional_data_handler=diffusion_data_handler,
target=target, streaming_pattern='pull', cpu_openmp=openmp, data_type=field_type_pdfs)
generate_boundary(ctx, 'thermal_LB_NeumannByCopy', NeumannByCopy(), method_thermal,
diff --git a/python/lbmpy_walberla/additional_data_handler.py b/python/lbmpy_walberla/additional_data_handler.py
index daaab32fd40aff5f2150495bf4782ef61a282dae..08208e39aba25a2813473bad49777c5b1e3c8e1e 100644
--- a/python/lbmpy_walberla/additional_data_handler.py
+++ b/python/lbmpy_walberla/additional_data_handler.py
@@ -9,7 +9,7 @@ try:
except ImportError:
from lbmpy.custom_code_nodes import MirroredStencilDirections
from lbmpy.boundaries.boundaryconditions import LbBoundary
-from lbmpy.boundaries import (ExtrapolationOutflow, FreeSlip, UBB, DiffusionDirichlet,
+from lbmpy.boundaries import (ExtrapolationOutflow, FreeSlip, UBB, DiffusionDirichlet, FixedDensity,
NoSlipLinearBouzidi, QuadraticBounceBack)
from pystencils_walberla.additional_data_handler import AdditionalDataHandler
@@ -41,6 +41,8 @@ def default_additional_data_handler(boundary_obj: LbBoundary, lb_method, field_n
return NoSlipLinearBouzidiAdditionalDataHandler(lb_method.stencil, boundary_obj)
elif isinstance(boundary_obj, QuadraticBounceBack):
return QuadraticBounceBackAdditionalDataHandler(lb_method.stencil, boundary_obj)
+ elif isinstance(boundary_obj, FixedDensity):
+ return FixedDensityAdditionalDataHandler(lb_method.stencil, boundary_obj)
else:
raise ValueError(f"No default AdditionalDataHandler available for boundary of type {boundary_obj.__class__}")
@@ -374,3 +376,38 @@ class DiffusionDirichletAdditionalDataHandler(AdditionalDataHandler):
def additional_member_variable(self):
return "std::function<real_t(const Cell &, const shared_ptr<StructuredBlockForest>&, IBlock&)> " \
"elementInitaliser; "
+
+
+class FixedDensityAdditionalDataHandler(AdditionalDataHandler):
+ def __init__(self, stencil, boundary_object):
+ assert isinstance(boundary_object, FixedDensity)
+ self._boundary_object = boundary_object
+ super(FixedDensityAdditionalDataHandler, self).__init__(stencil=stencil)
+
+ @property
+ def constructor_arguments(self):
+ return ", std::function<real_t(const Cell &, const shared_ptr<StructuredBlockForest>&, IBlock&)>& " \
+ "densityCallback "
+
+ @property
+ def initialiser_list(self):
+ return "elementInitaliser(densityCallback),"
+
+ @property
+ def additional_arguments_for_fill_function(self):
+ return "blocks, "
+
+ @property
+ def additional_parameters_for_fill_function(self):
+ return " const shared_ptr<StructuredBlockForest> &blocks, "
+
+ def data_initialisation(self, *_):
+ init_list = ["real_t InitialisatonAdditionalData = elementInitaliser(Cell(it.x(), it.y(), it.z()), "
+ "blocks, *block);", "element.rho = InitialisatonAdditionalData;"]
+
+ return "\n".join(init_list)
+
+ @property
+ def additional_member_variable(self):
+ return "std::function<real_t(const Cell &, const shared_ptr<StructuredBlockForest>&, IBlock&)> " \
+ "elementInitaliser; "