diff --git a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/3D_Validations.cpp b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/3D_Validations.cpp
index e784d166d41a33e2d63dde5aad7f763421d709e4..30bcad8f599c190062784f651464b62f060a5e12 100644
--- a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/3D_Validations.cpp
+++ b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/3D_Validations.cpp
@@ -161,7 +161,7 @@ int main(int argc, char** argv)
const real_t dt_SI = real_c((dx_SI * dx_SI * kappa_L) / (kappa_SI));
const real_t omegaConcentration = lbm::collision_model::omegaFromViscosity(kappa_L);
- Vector3< uint_t > domainSize(uint_c((xSize_SI / dx_SI)), uint_c((ySize_SI / dx_SI)), uint_c((zSize_SI/1)));
+ Vector3< uint_t > domainSize(uint_c((xSize_SI / dx_SI)), uint_c((ySize_SI / dx_SI)), uint_c((zSize_SI/dx_SI)));
const uint_t timeSteps = uint_c(std::ceil(tSI / dt_SI));
const uint_t tinitial = uint_c(0);
const real_t sigma_0 = real_t(10);
@@ -322,7 +322,7 @@ int main(int argc, char** argv)
// map boundaries into the fluid field simulation
lbm::BC_Fluid_Density density_fluid_bc(blocks, pdfFieldFluidCPUGPUID, real_t(1.0));
lbm::BC_Fluid_NoSlip noSlip_fluid_bc(blocks, pdfFieldFluidCPUGPUID);
- lbm::BC_Fluid_UBB ubb_fluid_bc(blocks, pdfFieldFluidCPUGPUID, uInflow_LBM, real_t(0));
+ lbm::BC_Fluid_UBB ubb_fluid_bc(blocks, pdfFieldFluidCPUGPUID, uInflow_LBM, real_t(0),real_t(0));
lbm::BC_Fluid_Neumann neumann_fluid_bc(blocks, pdfFieldFluidCPUGPUID);
lbm::BC_Concentration_Density density_concentration_bc_west(blocks, pdfFieldConcentrationCPUGPUID, real_t(0));
lbm::BC_Concentration_Density density_concentration_bc_south(blocks, pdfFieldConcentrationCPUGPUID, real_t(1));
@@ -345,12 +345,15 @@ int main(int argc, char** argv)
#else
- initConcentrationFieldGaussian(blocks, densityConcentrationFieldCPUGPUID, simulationDomain, domainSize, sigma_0,
- sigma_D, Uinitialize, x_0);
+ initConcentrationFieldSinusoidal(blocks,
+ densityConcentrationFieldCPUGPUID,simulationDomain,
+ domainSize,sigma_0,sigma_D,
+ Uinitialize,x_0,dx_SI,dt_SI);
// initConcentrationField(blocks,densityConcentrationFieldCPUGPUID,simulationDomain,domainSize);
-
- initFluidField(blocks, velFieldFluidCPUGPUID, Uinitialize);
+ WALBERLA_LOG_INFO_ON_ROOT("reached hereeee");
+ initFluidField(blocks, velFieldFluidCPUGPUID, Uinitialize,domainSize);
+ WALBERLA_LOG_INFO_ON_ROOT("reached hereeee 22");
#endif
@@ -516,18 +519,19 @@ int main(int argc, char** argv)
WcTimingPool timeloopTiming;
// TODO: maybe add warmup phase
- real_t advection_time = domainSize[0]/Uinitialize[0];
- analyticalSolGaussian(blocks, analyticalConcentrationFieldID, simulationDomain, domainSize, sigma_0, diffusivity,
- Uinitialize, x_0, 0,1);
+
+ analyticalSolSinusoidal(blocks,
+ analyticalConcentrationFieldID,simulationDomain,
+ domainSize, sigma_0,diffusivity,
+ Uinitialize,x_0,0,dx_SI,dt_SI);
for (uint_t timeStep = 1; timeStep <= timeSteps; ++timeStep)
{
if (useCommunicationHiding) { commTimeloop.singleStep(timeloopTiming); }
timeloop.singleStep(timeloopTiming);
-
- uint_t advection_period = uint_c(timeStep/advection_time);
- WALBERLA_LOG_INFO("advection period is " << timeStep << " " << advection_period);
- analyticalSolGaussian(blocks, analyticalConcentrationFieldID, simulationDomain, domainSize, sigma_0, diffusivity,
- Uinitialize, x_0, timeStep,advection_period);
+ analyticalSolSinusoidal(blocks,
+ analyticalConcentrationFieldID,simulationDomain,
+ domainSize, sigma_0,diffusivity,
+ Uinitialize,x_0,timeStep,dx_SI,dt_SI);
std::vector< real_t > norms = computeErrorL2(blocks, densityConcentrationFieldCPUGPUID,
analyticalConcentrationFieldID, ErrorFieldID, simulationDomain);
WALBERLA_LOG_INFO_ON_ROOT("Infinity norm is " << norms[0] << "\n L1 norm is " << norms[1] << "\n L2 norm is "
diff --git a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/CMakeLists.txt b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/CMakeLists.txt
index 1c5bb8aa4674af265990956ee8d2b6a082c7bd80..eb1db7410958fd0eea3c528bfa4cf9d1a9a270ea 100644
--- a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/CMakeLists.txt
+++ b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/CMakeLists.txt
@@ -33,8 +33,8 @@ if (WALBERLA_BUILD_WITH_CODEGEN)
GeneralInfoHeader.h
)
- waLBerla_add_executable(NAME 3d_gaussian_paper FILES 3D_Validations.cpp ../../utilities/InitializerFunctions.cpp
+ waLBerla_add_executable(NAME MTW_1way_3Dvalidations FILES 3D_Validations.cpp ../../utilities/InitializerFunctions.cpp
DEPENDS blockforest core field geometry gpu lbm lbm_mesapd_coupling mesa_pd sqlite vtk Validations_3D_${config})
- target_compile_definitions(3d_gaussian_paper PRIVATE Weighting=2)
+ target_compile_definitions(MTW_1way_3Dvalidations PRIVATE Weighting=2)
endif ()
endif ()
diff --git a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/Validations_3D.py b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/Validations_3D.py
index 021c768f8aa664adc7f5c4cd4e4b726489cad579..412f624118b2c9a9291c58490a992a0604480515 100644
--- a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/Validations_3D.py
+++ b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/Validations_3D.py
@@ -4,7 +4,7 @@ import pystencils as ps
from sympy.core.add import Add
from sympy.codegen.ast import Assignment
import sys
-from lbmpy.session import *
+#from lbmpy.session import *
from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil, ForceModel
from lbmpy.partially_saturated_cells import PSMConfig
@@ -41,13 +41,13 @@ const bool infoCsePdfs = {cse_pdfs};
with CodeGeneration() as ctx:
data_type = "float64" if ctx.double_accuracy else "float32"
- stencil_fluid = LBStencil(Stencil.D2Q9)
- stencil_concentration = LBStencil(Stencil.D2Q9)
+ stencil_fluid = LBStencil(Stencil.D3Q19)
+ stencil_concentration = LBStencil(Stencil.D3Q19)
omega = sp.Symbol("omega") # for now same for both the sweeps
omega_c = sp.Symbol("omega_c")
init_density_fluid = sp.Symbol("init_density_fluid")
init_density_concentration = sp.Symbol("init_density_concentration")
- init_velocity_fluid = sp.symbols("init_velocity_fluid_:2")
+ init_velocity_fluid = sp.symbols("init_velocity_fluid_:3")
#init_velocity_concentration = sp.symbols("init_velocity_concentration_:3")
pdfs_inter_fluid = sp.symbols("pdfs_inter_fluid:" + str(stencil_fluid.Q))
pdfs_inter_concentration = sp.symbols("pdfs_inter_concentration:" + str(stencil_concentration.Q))
@@ -66,13 +66,13 @@ with CodeGeneration() as ctx:
# Fluid PDFs and fields
pdfs_fluid, pdfs_fluid_tmp, velocity_field, density_field = ps.fields(
- f"pdfs_fluid({stencil_fluid.Q}), pdfs_fluid_tmp({stencil_fluid.Q}), velocity_field({stencil_fluid.D}), density_field({1}): {data_type}[2D]",
+ f"pdfs_fluid({stencil_fluid.Q}), pdfs_fluid_tmp({stencil_fluid.Q}), velocity_field({stencil_fluid.D}), density_field({1}): {data_type}[3D]",
layout=layout,
)
# Concentration PDFs and fields
pdfs_concentration, pdfs_concentration_tmp, concentration_field = ps.fields(
- f"pdfs_concentration({stencil_concentration.Q}), pdfs_concentration_tmp({stencil_concentration.Q}), concentration_field({1}): {data_type}[2D]",
+ f"pdfs_concentration({stencil_concentration.Q}), pdfs_concentration_tmp({stencil_concentration.Q}), concentration_field({1}): {data_type}[3D]",
layout=layout,
)
@@ -80,7 +80,7 @@ with CodeGeneration() as ctx:
if config_tokens[1]== "1":
concentration_output = None
- force_on_fluid = sp.symbols("F_:2")
+ force_on_fluid = sp.symbols("F_:3")
print("One-way fluid-concentration coupling set")
elif config_tokens[1] == "2":
@@ -122,7 +122,7 @@ with CodeGeneration() as ctx:
kernel_type="stream_pull_collide",
)
lamda=sp.Rational(1,6)
- omegacTRT = relaxation_rate_from_magic_number(omega_c, magic_number=lamda)
+ #omegacTRT = relaxation_rate_from_magic_number(omega_c, magic_number=lamda)
# Concentration LBM config
lbm_concentration_config = LBMConfig(
stencil=stencil_concentration,
@@ -227,7 +227,7 @@ with CodeGeneration() as ctx:
target=target,
)
- bc_velocity_fluid = sp.symbols("bc_velocity_fluid_:2")
+ bc_velocity_fluid = sp.symbols("bc_velocity_fluid_:3")
generate_boundary(
ctx,
"BC_Fluid_UBB",
@@ -281,7 +281,7 @@ with CodeGeneration() as ctx:
target=target,
)
- bc_velocity_concentration = sp.symbols("bc_velocity_concentration_:2") ## is it needed ?
+ bc_velocity_concentration = sp.symbols("bc_velocity_concentration_:3") ## is it needed ?
generate_boundary(
ctx,
"BC_Concentration_UBB",
diff --git a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/parameter_files/3d_gaussian.prm b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/parameter_files/3d_gaussian.prm
index 8578df8f18351cf4efc1197db9ef13ef28fa7f64..3630fbe03f5d333d0d45fe8acad8cdf5c49e5e63 100644
--- a/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/parameter_files/3d_gaussian.prm
+++ b/apps/benchmarks/MaterialTransport/one_way_coupling/3D_Validations/parameter_files/3d_gaussian.prm
@@ -1,13 +1,13 @@
Simulation{
- simName 3D_resolution; // 3D_validation or 3D_resolution
+ simName 3D_validation; // 3D_validation or 3D_resolution
}
NumericalSetup
{
- xSize 256; // SI and LBM length are one and the same => dx_SI = dx_LB always
- ySize 256;
- zSize 1;
+ xSize 2; // SI and LBM length are one and the same => dx_SI = dx_LB always
+ ySize 2;
+ zSize 2;
numXBlocks 1;
numYBlocks 1;
@@ -24,10 +24,10 @@ NumericalSetup
frameWidth <1, 1, 1>; // width of the outer region if splitting the LBM/PSM into inner and outer (only used if useCommunicationHiding is true)
- dx_SI 1;
+ dx_SI 0.02;
// fluid quantities
- uInflowL 0.8; // LBM velocity of fluid
+ uInflowL 1e-8; // LBM velocity of fluid
kappa_SI 4; // diffusivity SI
kappa_L 0.00000166; // diffusivity lbm
}
diff --git a/apps/benchmarks/MaterialTransport/utilities/InitializerFunctions.cpp b/apps/benchmarks/MaterialTransport/utilities/InitializerFunctions.cpp
index a77c5c0159049bb9ada2c0254d264fcae288d528..c210a8d42f011ab0ae4854aef1caf4d71f035e45 100644
--- a/apps/benchmarks/MaterialTransport/utilities/InitializerFunctions.cpp
+++ b/apps/benchmarks/MaterialTransport/utilities/InitializerFunctions.cpp
@@ -108,7 +108,7 @@ void initConcentrationFieldSinusoidal(const shared_ptr< StructuredBlockStorage >
blocks->transformBlockLocalToGlobalCell(globalCell, block, *cellIt);
Vector3< real_t > pos = blocks->getCellCenter(globalCell, level);
- ConcentrationField->get(*cellIt) = std::sin(pi*pos[0]*dx)*std::sin(pi*pos[1]*dx)*std::sin(pi*pos[2]*dx) +1;
+ ConcentrationField->get(*cellIt) = std::sin(pi*pos[0]*dx)*std::sin(pi*pos[1]*dx)*std::sin(pi*pos[2]*dx) +1 ;
ConcentrationField->get(*cellIt) = std::max(ConcentrationField->get(*cellIt), 1e-15);
}
}
@@ -187,7 +187,7 @@ void initFluidFieldPoiseuille(const shared_ptr< StructuredBlockStorage >& blocks
void analyticalSolGaussian(const shared_ptr< StructuredBlockStorage >& blocks,
BlockDataID& AnalyticalConcentrationFieldID, const math::AABB& domainAABB,
Vector3< uint_t > domainSize, const real_t sigma_0, const real_t diffusivity,
- const Vector3< real_t > uInflow, const Vector3< real_t > x_0, const real_t time)
+ const Vector3< real_t > uInflow, const Vector3< real_t > x_0, const real_t time,uint_t advection_period)
{
real_t sigma_D2 = (2 * diffusivity * time);