Skip to content

Commit

Permalink
add optional initial guess defaults for velocity and scalars
Browse files Browse the repository at this point in the history
  • Loading branch information
Stefan Kerkemeier committed Oct 31, 2020
1 parent 10b3fed commit 2c6f01f
Show file tree
Hide file tree
Showing 7 changed files with 81 additions and 24 deletions.
43 changes: 34 additions & 9 deletions okl/core/insVelocityExt.okl → okl/core/insExtrapolate.okl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,23 +24,23 @@

*/

@kernel void insVelocityExt(const dlong Nelements,
const int Nstages,
const dlong fieldOffset,
@restrict const dfloat* c,
@restrict const dfloat* U,
@restrict dfloat* Ue)
@kernel void insMultiExtrapolate(const dlong Nelements,
const int Nfields,
const int Nstages,
const dlong fieldOffset,
@restrict const dfloat* c,
@restrict const dfloat* U,
@restrict dfloat* Ue)
{
for(dlong e = 0; e < Nelements; ++e; @outer(0)) {
for(int n = 0; n < p_Np; ++n; @inner(0)) {
const dlong id = n + p_Np * e;

#pragma unroll p_NVfields
for (int i = 0; i < p_NVfields; i++) {
for (int i = 0; i < Nfields; i++) {
dfloat Un = 0.;

for (int s = 0; s < Nstages; s++) {
const dlong idm = id + i * fieldOffset + s * p_NVfields * fieldOffset;
const dlong idm = id + i * fieldOffset + s * Nfields * fieldOffset;
const dfloat Um = U[idm];
Un += c[s] * Um;
}
Expand All @@ -50,3 +50,28 @@
}
}
}

@kernel void insExtrapolate(const dlong Nelements,
const int Nstages,
const dlong fieldOffset,
const dlong stageOffset,
@restrict const dfloat* c,
@restrict const dfloat* U,
@restrict dfloat* Ue)
{
for(dlong e = 0; e < Nelements; ++e; @outer(0)) {
for(int n = 0; n < p_Np; ++n; @inner(0)) {
const dlong id = n + p_Np * e;

dfloat Un = 0.;

for (int s = 0; s < Nstages; s++) {
const dlong idm = id + fieldOffset + s * stageOffset;
const dfloat Um = U[idm];
Un += c[s] * Um;
}

Ue[id] = Un;
}
}
}
6 changes: 5 additions & 1 deletion src/core/cds.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@ occa::memory cdsSolve(const int is, cds_t* cds, dfloat time)
}
elliptic_t* solver = cds->solver[is];

//copy current solution fields as initial guess
cds->o_wrk0.copyFrom(cds->o_S, cds->Ntotal * sizeof(dfloat), 0, is * cds->fieldOffset * sizeof(dfloat));

//enforce Dirichlet BCs
Expand Down Expand Up @@ -60,6 +59,11 @@ occa::memory cdsSolve(const int is, cds_t* cds, dfloat time)
oogs::startFinish(cds->o_wrk1, 1, cds->fieldOffset, ogsDfloat, ogsAdd, gsh);
if (solver->Nmasked) mesh->maskKernel(solver->Nmasked, solver->o_maskIds, cds->o_wrk1);

if(cds->options.compareArgs("SCALAR INITIAL GUESS DEFAULT", "EXTRAPOLATION")) {
cds->o_wrk0.copyFrom(cds->o_Se, cds->Ntotal * sizeof(dfloat), 0, is * cds->fieldOffset * sizeof(dfloat));
if (solver->Nmasked) cds->maskCopyKernel(solver->Nmasked, 0, solver->o_maskIds, cds->o_wrk2, cds->o_wrk0);
}

cds->Niter[is] = ellipticSolve(solver, cds->TOL, cds->o_wrk1, cds->o_wrk0);

return cds->o_wrk0;
Expand Down
3 changes: 2 additions & 1 deletion src/core/cds.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,6 +105,7 @@ typedef struct

occa::memory o_wrk0, o_wrk1, o_wrk2, o_wrk3, o_wrk4, o_wrk5, o_wrk6;

occa::kernel extrapolateKernel;
occa::kernel fillKernel;
occa::kernel sumMakefKernel;
occa::kernel scaledAddKernel;
Expand All @@ -119,7 +120,7 @@ typedef struct
// occa::kernel constrainKernel;

occa::memory o_U;
occa::memory o_S;
occa::memory o_S, o_Se;

// occa::memory o_Vort, o_Div; // Not sure to keep it
occa::memory o_haloBuffer;
Expand Down
4 changes: 1 addition & 3 deletions src/core/ins.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,6 @@ typedef struct

dfloat* U, * P;
dfloat* BF, * FU;
dfloat* PI;

//RK Subcycle Data
int SNrk;
Expand Down Expand Up @@ -122,7 +121,7 @@ typedef struct
occa::kernel subCycleVolumeKernel, subCycleCubatureVolumeKernel;
occa::kernel subCycleSurfaceKernel, subCycleCubatureSurfaceKernel;
occa::kernel subCycleRKUpdateKernel;
occa::kernel velocityExtKernel;
occa::kernel extrapolateKernel;

occa::kernel wgradientVolumeKernel;

Expand All @@ -141,7 +140,6 @@ typedef struct
occa::memory o_prop, o_ellipticCoeff;

occa::memory o_UH;
occa::memory o_PI;

occa::memory o_vHaloBuffer, o_pHaloBuffer;
occa::memory o_velocityHaloGatherTmp;
Expand Down
18 changes: 11 additions & 7 deletions src/core/insSetup.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,6 @@ ins_t* insSetup(MPI_Comm comm, occa::device device, setupAide &options, int buil
ins->Ue = (dfloat*) calloc(ins->NVfields * ins->fieldOffset,sizeof(dfloat));

ins->P = (dfloat*) calloc(ins->fieldOffset,sizeof(dfloat));
ins->PI = (dfloat*) calloc(ins->fieldOffset,sizeof(dfloat));

ins->BF = (dfloat*) calloc(ins->NVfields * ins->fieldOffset,sizeof(dfloat));
ins->FU = (dfloat*) calloc(ins->NVfields * ins->Nstages * ins->fieldOffset,sizeof(dfloat));
Expand Down Expand Up @@ -189,7 +188,6 @@ ins_t* insSetup(MPI_Comm comm, occa::device device, setupAide &options, int buil
ins->U);
ins->o_Ue = mesh->device.malloc(ins->NVfields * ins->fieldOffset * sizeof(dfloat), ins->Ue);
ins->o_P = mesh->device.malloc(ins->fieldOffset * sizeof(dfloat), ins->P);
ins->o_PI = mesh->device.malloc(ins->fieldOffset * sizeof(dfloat), ins->PI);

ins->o_FU =
mesh->device.malloc(ins->NVfields * ins->Nstages * ins->fieldOffset * sizeof(dfloat),
Expand Down Expand Up @@ -712,9 +710,9 @@ ins_t* insSetup(MPI_Comm comm, occa::device device, setupAide &options, int buil
ins->subCycleRKUpdateKernel =
mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);

fileName = oklpath + "insVelocityExt" + ".okl";
kernelName = "insVelocityExt";
ins->velocityExtKernel =
fileName = oklpath + "insExtrapolate" + ".okl";
kernelName = "insMultiExtrapolate";
ins->extrapolateKernel =
mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);

// ===========================================================================
Expand Down Expand Up @@ -947,8 +945,9 @@ cds_t* cdsSetup(ins_t* ins, mesh_t* mesh, setupAide options, occa::properties &k
cds->o_U = ins->o_U;
cds->o_Ue = ins->o_Ue;
cds->o_S =
mesh->device.malloc(cds->NSfields * cds->Nstages * cds->fieldOffset * sizeof(dfloat),
cds->S);
mesh->device.malloc(cds->NSfields * cds->Nstages * cds->fieldOffset * sizeof(dfloat), cds->S);
cds->o_Se =
mesh->device.malloc(cds->NSfields * cds->Nstages * cds->fieldOffset * sizeof(dfloat), cds->S);
cds->o_BF = mesh->device.malloc(cds->NSfields * cds->fieldOffset * sizeof(dfloat), cds->BF);
cds->o_FS =
mesh->device.malloc(cds->NSfields * cds->Nstages * cds->fieldOffset * sizeof(dfloat),
Expand Down Expand Up @@ -1145,6 +1144,11 @@ cds_t* cdsSetup(ins_t* ins, mesh_t* mesh, setupAide options, occa::properties &k
cds->scaledAddKernel =
mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);

fileName = oklpath + "insExtrapolate" + ".okl";
kernelName = "insExtrapolate";
cds->extrapolateKernel =
mesh->device.buildKernel(fileName.c_str(), kernelName.c_str(), kernelInfo);

if(cds->Nsubsteps) {
fileName = oklpath + "cdsSubCycle" + suffix + ".okl";
kernelName = "cdsSubCycleStrongCubatureVolume" + suffix;
Expand Down
13 changes: 11 additions & 2 deletions src/core/runTime.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,14 +41,23 @@ void runStep(ins_t* ins, dfloat time, dfloat dt, int tstep)
else if(tstep <= 3 && ins->temporalOrder >= 3)
extbdfCoefficents(ins,tstep);

// First extrapolate velocity to t^(n+1)
// extrapolate
if(ins->flow)
ins->velocityExtKernel(mesh->Nelements,
ins->extrapolateKernel(mesh->Nelements,
ins->NVfields,
ins->ExplicitOrder,
ins->fieldOffset,
ins->o_extbdfA,
ins->o_U,
ins->o_Ue);
if(ins->Nscalar)
ins->extrapolateKernel(mesh->Nelements,
cds->NSfields,
cds->ExplicitOrder,
cds->fieldOffset,
cds->o_extbdfA,
cds->o_S,
cds->o_Se);

if(ins->Nscalar)
scalarSolve(ins, time, dt, cds->o_S);
Expand Down
18 changes: 17 additions & 1 deletion src/core/tombo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,7 +249,23 @@ occa::memory velocitySolve(ins_t* ins, dfloat time)

oogs::startFinish(ins->o_wrk3, ins->NVfields, ins->fieldOffset,ogsDfloat, ogsAdd, ins->gsh);

ins->o_wrk0.copyFrom(ins->o_U, ins->NVfields * ins->fieldOffset * sizeof(dfloat));
if(ins->options.compareArgs("VELOCITY INITIAL GUESS DEFAULT", "EXTRAPOLATION")) {
ins->o_wrk0.copyFrom(ins->o_Ue, ins->NVfields * ins->fieldOffset * sizeof(dfloat));
if (ins->uvwSolver) {
if (ins->uvwSolver->Nmasked) ins->maskCopyKernel(ins->uvwSolver->Nmasked, 0*ins->fieldOffset, ins->uvwSolver->o_maskIds,
ins->o_U, ins->o_wrk0);
} else {
if (ins->uSolver->Nmasked) ins->maskCopyKernel(ins->uSolver->Nmasked, 0*ins->fieldOffset, ins->uSolver->o_maskIds,
ins->o_U, ins->o_wrk0);
if (ins->vSolver->Nmasked) ins->maskCopyKernel(ins->vSolver->Nmasked, 1*ins->fieldOffset, ins->vSolver->o_maskIds,
ins->o_U, ins->o_wrk0);
if (ins->wSolver->Nmasked) ins->maskCopyKernel(ins->wSolver->Nmasked, 2*ins->fieldOffset, ins->wSolver->o_maskIds,
ins->o_U, ins->o_wrk0);
}
} else {
ins->o_wrk0.copyFrom(ins->o_U, ins->NVfields * ins->fieldOffset * sizeof(dfloat));
}

if(ins->uvwSolver) {
ins->NiterU = ellipticSolve(ins->uvwSolver, ins->velTOL, ins->o_wrk3, ins->o_wrk0);
} else {
Expand Down

0 comments on commit 2c6f01f

Please sign in to comment.