From 2c6f01f612f3f358bbc62d0ee3442450d3ca1c79 Mon Sep 17 00:00:00 2001 From: Stefan Kerkemeier Date: Sat, 31 Oct 2020 08:47:19 +0000 Subject: [PATCH] add optional initial guess defaults for velocity and scalars --- ...{insVelocityExt.okl => insExtrapolate.okl} | 43 +++++++++++++++---- src/core/cds.cpp | 6 ++- src/core/cds.h | 3 +- src/core/ins.h | 4 +- src/core/insSetup.cpp | 18 +++++--- src/core/runTime.cpp | 13 +++++- src/core/tombo.cpp | 18 +++++++- 7 files changed, 81 insertions(+), 24 deletions(-) rename okl/core/{insVelocityExt.okl => insExtrapolate.okl} (64%) diff --git a/okl/core/insVelocityExt.okl b/okl/core/insExtrapolate.okl similarity index 64% rename from okl/core/insVelocityExt.okl rename to okl/core/insExtrapolate.okl index add64c2a7..21bf1118a 100644 --- a/okl/core/insVelocityExt.okl +++ b/okl/core/insExtrapolate.okl @@ -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; } @@ -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; + } + } +} diff --git a/src/core/cds.cpp b/src/core/cds.cpp index 6ce841ee9..6e1a23a94 100644 --- a/src/core/cds.cpp +++ b/src/core/cds.cpp @@ -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 @@ -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; diff --git a/src/core/cds.h b/src/core/cds.h index a0d419c17..7e92f5644 100644 --- a/src/core/cds.h +++ b/src/core/cds.h @@ -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; @@ -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; diff --git a/src/core/ins.h b/src/core/ins.h index 718fa4080..32c8b84ba 100644 --- a/src/core/ins.h +++ b/src/core/ins.h @@ -63,7 +63,6 @@ typedef struct dfloat* U, * P; dfloat* BF, * FU; - dfloat* PI; //RK Subcycle Data int SNrk; @@ -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; @@ -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; diff --git a/src/core/insSetup.cpp b/src/core/insSetup.cpp index f4165c9ed..e1f730578 100644 --- a/src/core/insSetup.cpp +++ b/src/core/insSetup.cpp @@ -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)); @@ -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), @@ -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); // =========================================================================== @@ -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), @@ -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; diff --git a/src/core/runTime.cpp b/src/core/runTime.cpp index 8e8e5327c..d0f0b5504 100644 --- a/src/core/runTime.cpp +++ b/src/core/runTime.cpp @@ -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); diff --git a/src/core/tombo.cpp b/src/core/tombo.cpp index 2eadeb89a..ccba2907b 100644 --- a/src/core/tombo.cpp +++ b/src/core/tombo.cpp @@ -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 {