diff --git a/apps/gpu/RotatingGrid/RotatingGrid.cpp b/apps/gpu/RotatingGrid/RotatingGrid.cpp
index 0826fb2f1abcc4958ad9b269d712a58b689531ac..99c58f2e4c85bf1300f358a0181cfa15a1cc3dcc 100644
--- a/apps/gpu/RotatingGrid/RotatingGrid.cpp
+++ b/apps/gpu/RotatingGrid/RotatingGrid.cpp
@@ -91,8 +91,9 @@ int main()
const real velocityLB = 0.05; // LB units
const uint nx = 64;
- const uint timeStepOut = 1;
- const uint timeStepEnd = 10;
+ const uint timeStepOut = 1000;
+ const uint timeStepEnd = 10000;
+ const uint timeStepStartOutput = 0;
//////////////////////////////////////////////////////////////////////////
// compute parameters in lattice units
@@ -108,10 +109,10 @@ int main()
//////////////////////////////////////////////////////////////////////////
auto gridBuilder = std::make_shared<MultipleGridBuilder>();
- gridBuilder->addCoarseGrid(-1.5 * L, -0.5 * L, -0.5 * L, 1.5 * L, 0.5 * L, 0.5 * L, dx);
+ gridBuilder->addCoarseGrid(-0.5 * L, -0.5 * L, -0.5 * L, 2.0 * L, 0.5 * L, 0.5 * L, dx);
if (rotOrInt == Rot) gridBuilder->addGridRotatingGrid(std::make_shared<Cylinder>(0.1, 0.1, 0.1, 0.25 * L, 1. * L, Axis::x));
- if (rotOrInt == Int) gridBuilder->addGrid(std::make_shared<Cylinder>(0.0, 0.0, 0.0, 0.3 * L, 1. * L, Axis::x), 1);
+ if (rotOrInt == Int) gridBuilder->addGrid(std::make_shared<Cylinder>(0.1, 0.1, 0.1, 0.25 * L, 0.8 * L, Axis::x), 1);
GridScalingFactory scalingFactory = GridScalingFactory();
scalingFactory.setScalingFactory(GridScalingFactory::GridScaling::ScaleCompressible);
@@ -138,6 +139,7 @@ int main()
para->setTimestepOut(timeStepOut);
para->setTimestepEnd(timeStepEnd);
+ para->setTimestepStartOut(timeStepStartOutput);
para->configureMainKernel(vf::CollisionKernel::Compressible::K17CompressibleNavierStokes);
@@ -145,11 +147,18 @@ int main()
// set boundary conditions
//////////////////////////////////////////////////////////////////////////
- gridBuilder->setSlipBoundaryCondition(SideType::MY, 0.0, 0.0, 0.0);
- gridBuilder->setSlipBoundaryCondition(SideType::PY, 0.0, 0.0, 0.0);
- gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
- gridBuilder->setSlipBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
+ // gridBuilder->setSlipBoundaryCondition(SideType::MY, 0.0, 0.0, 0.0);
+ // gridBuilder->setSlipBoundaryCondition(SideType::PY, 0.0, 0.0, 0.0);
+ // gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
+ // gridBuilder->setSlipBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
+ gridBuilder->setNoSlipBoundaryCondition(SideType::MY);
+ gridBuilder->setNoSlipBoundaryCondition(SideType::PY);
+ gridBuilder->setNoSlipBoundaryCondition(SideType::MZ);
+ gridBuilder->setNoSlipBoundaryCondition(SideType::PZ);
+
+
+ // gridBuilder->setVelocityBoundaryCondition(SideType::MX, 0.0, 0.0, 0.0);
gridBuilder->setVelocityBoundaryCondition(SideType::MX, velocityLB, 0.0, 0.0);
gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
@@ -164,12 +173,14 @@ int main()
// set initial condition
//////////////////////////////////////////////////////////////////////////
- para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
- rho = (real)0.0;
- vx = velocityLB;
- vy = 0.0;
- vz = 0.0;
- });
+ // para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
+ // rho = (real)0.0;
+ // // if (coordX > -0.52 && coordY > 0.27 && coordZ > -0.1 && coordX < -0.49 && coordY < 0.3 && coordZ < 0.11) rho = 1e-5;
+ // // if (coordX > -0.34 && coordY > 0.27 && coordZ > -0.1 && coordX < -0.32 && coordY < 0.3 && coordZ < 0.11) rho = 1e-5;
+ // vx = 0.0;
+ // vy = 0.0;
+ // vz = 0.0;
+ // });
//////////////////////////////////////////////////////////////////////////
// set copy mesh to simulation
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
index c1359069dc92d5332471d0c5d9ee2547004322fb..0ac72a9c0d873bd24e91252b8812a72f402b92c8 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
+++ b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
@@ -63,7 +63,6 @@ void UpdateGrid27::updateGrid(int level, unsigned int t)
paraRot->parameterRotDevice->gridAngle[0] += paraRot->parameterRotDevice->angularVelocity[0];
paraRot->parameterRotDevice->gridAngle[1] += paraRot->parameterRotDevice->angularVelocity[1];
paraRot->parameterRotDevice->gridAngle[2] += paraRot->parameterRotDevice->angularVelocity[2];
- VF_LOG_DEBUG("gridAngleX = {}", paraRot->parameterRotDevice->gridAngle[0]);
rotationInterpolation(level);
} else {
refinement(this, para.get(), level);
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
index 6851e795c4453734a0bc12380bf0283d5942c60a..43489caa3132f30c76ed64d0dce089f146147ed5 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
@@ -119,6 +119,18 @@ void GridGenerator::allocArrays_CoordNeighborGeo()
cudaMemoryManager->cudaAllocCoordRotation(1);
para->fillCoordinateVectorsForRotatingGrid(1);
cudaMemoryManager->cudaCopyCoordRotation(1);
+ auto paraStaticGridH = para->getParH(0);
+ auto paraStaticGridD = para->getParD(0);
+ auto paraRotatingGridH = para->getParH(1);
+ auto paraRotatingGridD = para->getParD(1);
+ paraRotatingGridH->viscosity = paraStaticGridH->viscosity;
+ paraRotatingGridD->viscosity = paraStaticGridD->viscosity;
+ paraRotatingGridH->diffusivity = paraStaticGridH->diffusivity;
+ paraRotatingGridD->diffusivity = paraStaticGridD->diffusivity;
+ paraRotatingGridH->omega = paraStaticGridH->omega;
+ paraRotatingGridD->omega = paraStaticGridD->omega;
+ paraRotatingGridH->scalingFactorOfGridSpacing = 1.0;
+ paraRotatingGridD->scalingFactorOfGridSpacing = 1.0;
}
for (int i = 0; i <= para->getMaxLevel(); i++) {
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
index 4dbbb16194053ba019760378e093cc9607c761f4..bd870c83752168a492c77ddad4af3229fd48b1f5 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
@@ -174,8 +174,10 @@ __global__ void interpolateRotatingToStatic(
bool isEvenTimestep,
real dx)
{
+
const unsigned listIndex = vf::gpu::getNodeIndex();
if (listIndex >= numberOfInterfaceNodes) return;
+
const uint destinationIndex = indicesStatic[listIndex];
const uint previousSourceIndex = indicesRotatingCell[listIndex];
const uint indexNeighborMMMsource = neighborMMMrotating[previousSourceIndex];
@@ -201,7 +203,7 @@ __global__ void interpolateRotatingToStatic(
real forcesVertexX [8];
real forcesVertexY [8];
real forcesVertexZ [8];
- real tangentialVelocitiesX [8];
+ real tangentialVelocitiesX [8];
real tangentialVelocitiesY [8];
real tangentialVelocitiesZ [8];
@@ -274,16 +276,16 @@ __global__ void interpolateRotatingToStatic(
// 1. calculate moments for the nodes of the source cell, add half of the rotational force to the velocity during the computation of the moments
vf::lbm::MomentsOnSourceNodeSet momentsSet;
vf::gpu::calculateMomentSet<false>(momentsSet, listIndex, distributionsRotating, neighborXrotating, neighborYrotating,
- neighborZrotating, indicesRotatingCell, nullptr, numberOfLBNodesStatic, omegaRotating,
+ neighborZrotating, indicesRotatingCell, nullptr, numberOfLBNodesRotating, omegaRotating,
isEvenTimestep, forcesVertexX, forcesVertexY, forcesVertexZ);
// 2. calculate the coefficients for the interpolation
- // For this the relative coordinates of the destination point inside the source cell are needed
- // this cell coordinate system is centered at the middle of the source cell. The source cell has a side lenght of one
- // add tangential velocity to velocity for coefficient computaion
- const real cellCoordDestinationX = (coordSourceX[sourceIndex] + rotatedCoordDestinationX) / dx - c1o2;
- const real cellCoordDestinationY = (coordSourceY[sourceIndex] + rotatedCoordDestinationY) / dx - c1o2;
- const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] + rotatedCoordDestinationZ) / dx - c1o2;
+ // For this the relative coordinates of the source cell in the coordinate system of the destination node ("offsets")
+ // this cell coordinate system is centered at the destination node. The source cell has a side lenght of one.
+ // Also add the tangential velocity to velocity for coefficient computation
+ const real cellCoordDestinationX = (coordSourceX[sourceIndex] - rotatedCoordDestinationX) / dx + c1o2;
+ const real cellCoordDestinationY = (coordSourceY[sourceIndex] - rotatedCoordDestinationY) / dx + c1o2;
+ const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] - rotatedCoordDestinationZ) / dx + c1o2;
momentsSet.addToVelocity(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ);
vf::lbm::Coefficients coefficients;
momentsSet.calculateCoefficients(coefficients, cellCoordDestinationX, cellCoordDestinationY, cellCoordDestinationZ);
@@ -363,18 +365,21 @@ __global__ void interpolateRotatingToStatic(
vvyTemp = vvy;
vvzTemp = vvz;
- if (angleX != 0) {
+ if (angleX != c0o1) {
vvyTemp = vvy * cos(angleX) + vvz * sin(angleX);
vvzTemp = -vvy * sin(angleX) + vvz * cos(angleX);
- } else if (angleY != 0) {
+ } else if (angleY != c0o1) {
// rotate in y
vvxTemp = vvx * cos(angleY) - vvz * sin(angleY);
vvzTemp = vvx * sin(angleY) + vvz * cos(angleY);
- } else if (angleZ != 0) {
+ } else if (angleZ != c0o1) {
// rotate in z
vvxTemp = vvx * cos(angleZ) + vvy * sin(angleZ);
vvyTemp = -vvx * sin(angleZ) + vvy * cos(angleZ);
}
+ vvx = vvxTemp;
+ vvy = vvyTemp;
+ vvz = vvzTemp;
////////////////////////////////////////////////////////////////////////////////
// calculate the squares of the velocities
@@ -404,7 +409,7 @@ __global__ void interpolateRotatingToStatic(
m011 = m011SFOR;
m101 = m101SFOR;
m110 = m110SFOR;
- if (angleX != 0) {
+ if (angleX != c0o1) {
mxxMyy = (mxxMyySFOR + mxxMzzSFOR + (mxxMyySFOR - mxxMzzSFOR) * cos(angleX * c2o1) -
c2o1 * m011SFOR * sin(angleX * c2o1)) *
c1o2;
@@ -415,14 +420,14 @@ __global__ void interpolateRotatingToStatic(
m011 = m011SFOR * cos(angleX * c2o1) + (mxxMyySFOR - mxxMzzSFOR) * cos(angleX) * sin(angleX);
m101 = m101SFOR * cos(angleX) - m110SFOR * sin(angleX);
m110 = m110SFOR * cos(angleX) + m101SFOR * sin(angleX);
- } else if (angleY != 0) {
+ } else if (angleY != c0o1) {
mxxMyy = mxxMyySFOR - mxxMzzSFOR * c1o2 + c1o2 * mxxMzzSFOR * cos(angleY * c2o1) + m101SFOR * sin(angleY * c2o1);
mxxMzz = mxxMzzSFOR * cos(angleY * c2o1) + c2o1 * m101SFOR * sin(angleY * c2o1);
m011 = m011SFOR * cos(angleY) - m110SFOR * sin(angleY);
m101 = m101SFOR * cos(angleY * c2o1) - mxxMzzSFOR * cos(angleY) * sin(angleY);
m110 = m110SFOR * cos(angleY) + m011SFOR * sin(angleY);
- } else if (angleZ != 0) {
+ } else if (angleZ != c0o1) {
mxxMyy = mxxMyySFOR * cos(angleY * c2o1) - c2o1 * m110SFOR * sin(angleY * c2o1);
mxxMzz = -c1o2 * mxxMyySFOR + mxxMzzSFOR + c1o2 * mxxMyySFOR * cos(angleY * c2o1) - m110SFOR * sin(angleY * c2o1);
@@ -545,8 +550,8 @@ __global__ void interpolateRotatingToStatic(
vf::gpu::getPointersToDistributions(distStatic, distributionsStatic, numberOfLBNodesStatic, isEvenTimestep);
// write
- vf::gpu::ListIndices writeIndicesStatic(destinationIndex, neighborXstatic, neighborYstatic, neighborZstatic);
- vf::gpu::write(distStatic, writeIndicesStatic, fStatic);
+ vf::gpu::ListIndices indicesStaticForWriting(destinationIndex, neighborXstatic, neighborYstatic, neighborZstatic);
+ vf::gpu::write(distStatic, indicesStaticForWriting, fStatic);
}
__global__ void updateGlobalCoordinates(
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
index 70dd6ff823bbe6d56af8172ff1d544aceb4f7590..7707a1515ab84788701ad63e76b94d76520ae5ca 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
@@ -111,11 +111,11 @@ __global__ void interpolateStaticToRotating(
isEvenTimestep);
// 2. calculate the coefficients for the interpolation
- // For this the relative coordinates of the destination pooint inside the source cell are needed
- // this cell coordinate system is centered at the middle of the source cell. The source cell has a side lenght of one
- const real cellCoordDestinationX = (coordSourceX[sourceIndex] + globalCoordDestinationX) / dx - c1o2;
- const real cellCoordDestinationY = (coordSourceY[sourceIndex] + globalCoordDestinationY) / dx - c1o2;
- const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] + globalCoordDestinationZ) / dx - c1o2;
+ // For this the relative coordinates of the source cell in the coordinate system of the destination node ("offsets")
+ // this cell coordinate system is centered at the destination node. The source cell has a side lenght of one.
+ const real cellCoordDestinationX = (coordSourceX[sourceIndex] - globalCoordDestinationX) / dx + c1o2;
+ const real cellCoordDestinationY = (coordSourceY[sourceIndex] - globalCoordDestinationY) / dx + c1o2;
+ const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] - globalCoordDestinationZ) / dx + c1o2;
vf::lbm::Coefficients coefficients;
momentsSet.calculateCoefficients(coefficients, cellCoordDestinationX, cellCoordDestinationY, cellCoordDestinationZ);
@@ -185,7 +185,7 @@ __global__ void interpolateStaticToRotating(
//! - Declare local variables for destination nodes
//!
real vvx, vvy, vvz, vxsq, vysq, vzsq;
- real useNEQ = c1o1; // zero; //one; //.... one = on ..... zero = off
+ real useNEQ = c1o1; // c0o1; // c1o1; //.... one = on ..... zero = off
////////////////////////////////////////////////////////////////////////////////
//! - Set macroscopic values on destination node (zeroth and first order moments)
@@ -224,9 +224,9 @@ __global__ void interpolateStaticToRotating(
angularVelocityY * angularVelocityY * coordDestinationZlocal;
// Coriolis force
- forceX += c2o1 * (angularVelocityZ * vvy - angularVelocityY * vvz);
+ forceX += c2o1 * (angularVelocityZ * vvy - angularVelocityY * vvz);
forceY += -c2o1 * (angularVelocityZ * vvx - angularVelocityX * vvz);
- forceZ += c2o1 * (angularVelocityY * vvx - angularVelocityX * vvy);
+ forceZ += c2o1 * (angularVelocityY * vvx - angularVelocityX * vvy);
////////////////////////////////////////////////////////////////////////////////
//! - subtract half the force from the velocities
@@ -239,15 +239,17 @@ __global__ void interpolateStaticToRotating(
//! - rotate the velocities
//!
- real vvxTemp, vvyTemp, vvzTemp;
- if (angleX != 0) {
+ real vvxTemp = vvx;
+ real vvyTemp = vvy;
+ real vvzTemp = vvz;
+ if (angleX != c0o1) {
vvyTemp = vvy * cos(angleX) - vvz * sin(angleX);
vvzTemp = vvy * sin(angleX) + vvz * cos(angleX);
- } else if (angleY != 0) {
+ } else if (angleY != c0o1) {
// rotate in y
vvxTemp = vvx * cos(angleY) + vvz * sin(angleY);
vvzTemp = -vvx * sin(angleY) + vvz * cos(angleY);
- } else if (angleZ != 0) {
+ } else if (angleZ != c0o1) {
// rotate in z
vvxTemp = vvx * cos(angleZ) - vvy * sin(angleZ);
vvyTemp = vvx * sin(angleZ) + vvy * cos(angleZ);
@@ -284,7 +286,7 @@ __global__ void interpolateStaticToRotating(
m011 = m011SFOR;
m101 = m101SFOR;
m110 = m110SFOR;
- if (angleX != 0) {
+ if (angleX != c0o1) {
mxxMyy = c1o2 * (mxxMyySFOR + mxxMzzSFOR + ( mxxMyySFOR - mxxMzzSFOR) * cos(angleX * c2o1) +
c2o1 * m011SFOR * sin(angleX * c2o1));
mxxMzz = c1o2 * (mxxMyySFOR + mxxMzzSFOR + (-mxxMyySFOR + mxxMzzSFOR) * cos(angleX * c2o1) -
@@ -293,14 +295,14 @@ __global__ void interpolateStaticToRotating(
m011 = m011SFOR * cos(angleX * c2o1) + (-mxxMyySFOR + mxxMzzSFOR) * cos(angleX) * sin(angleX);
m101 = m101SFOR * cos(angleX) + m110SFOR * sin(angleX);
m110 = m110SFOR * cos(angleX) - m101SFOR * sin(angleX);
- } else if (angleY != 0) {
+ } else if (angleY != c0o1) {
mxxMyy = mxxMyySFOR - c1o2 * mxxMzzSFOR + c1o2 * mxxMzzSFOR * cos(angleY * c2o1) + m101SFOR * sin(angleY * c2o1);
mxxMzz = mxxMzzSFOR * cos(angleY * c2o1) + c2o1 * m101SFOR * sin(angleY * c2o1);
m011 = m011SFOR * cos(angleY) - m110SFOR * sin(angleY);
m101 = m101SFOR * cos(angleY * c2o1) - mxxMzzSFOR * cos(angleY) * sin(angleY);
m110 = m110SFOR * cos(angleY) + m011SFOR * sin(angleY);
- } else if (angleZ != 0) {
+ } else if (angleZ != c0o1) {
mxxMyy = mxxMyySFOR * cos(angleY * c2o1) - c2o1 * m110SFOR * sin(angleY * c2o1);
mxxMzz = -c1o2 * mxxMyySFOR + mxxMzzSFOR + c1o2 * mxxMyySFOR * cos(angleY * c2o1) - m110SFOR * sin(angleY * c2o1);
@@ -309,7 +311,7 @@ __global__ void interpolateStaticToRotating(
m110 = m110SFOR * cos(angleZ * c2o1) + mxxMyySFOR * cos(angleZ) * sin(angleZ);
}
- // calculate remaining second order moments from previously rotated moments
+ // calculate the remaining second order moments from previously rotated moments
m200 = c1o3 * ( mxxMyy + mxxMzz + mxxPyyPzz) * useNEQ;
m020 = c1o3 * (-c2o1 * mxxMyy + mxxMzz + mxxPyyPzz) * useNEQ;
m002 = c1o3 * ( mxxMyy - c2o1 * mxxMzz + mxxPyyPzz) * useNEQ;
@@ -422,6 +424,6 @@ __global__ void interpolateStaticToRotating(
vf::gpu::getPointersToDistributions(distRoating, distributionsRotating, numberOfLBNodesRotating, isEvenTimestep);
// write
- vf::gpu::ListIndices writeIndicesRotating(destinationIndex, neighborXrotating, neighborYrotating, neighborZrotating);
- vf::gpu::write(distRoating, writeIndicesRotating, fRotating);
+ vf::gpu::ListIndices indicesRotatingForWriting(destinationIndex, neighborXrotating, neighborYrotating, neighborZrotating);
+ vf::gpu::write(distRoating, indicesRotatingForWriting, fRotating);
}
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index ac5cd5c6d2079a4c19adb082cd182ead10a526ad..948f923c183c561f9b07de79b1d35308c70a00f1 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -972,10 +972,10 @@ void Simulation::readAndWriteFiles(uint timestep)
UpdateGlobalCoordinates(para->getParD(1).get(), para->getRotatingGridParameter()->parameterRotDevice.get());
cudaMemoryManager->cudaCopyCoordDeviceToHost(1);
// cudaMemoryManager->cudaCopyCoordRotationDeviceToHost(1);
- cudaMemoryManager->cudaCopyInterfaceCFDeviceToHost(0);
- cudaMemoryManager->cudaCopyInterfaceFCDeviceToHost(0);
- InterfaceDebugWriter::writeInterfaceLinesDebugCF(para.get(), timestep);
- InterfaceDebugWriter::writeInterfaceLinesDebugFC(para.get(), timestep);
+ // cudaMemoryManager->cudaCopyInterfaceCFDeviceToHost(0);
+ // cudaMemoryManager->cudaCopyInterfaceFCDeviceToHost(0);
+ // InterfaceDebugWriter::writeInterfaceLinesDebugCF(para.get(), timestep);
+ // InterfaceDebugWriter::writeInterfaceLinesDebugFC(para.get(), timestep);
////////////////////////////////////////////////////////////////////////