Commit aa9467e4 authored by Stefan Kerkemeier's avatar Stefan Kerkemeier
Browse files

compute residual same as in nek

parent 8a798e2e
......@@ -263,6 +263,7 @@
@restrict const dlong* vmapM,
@restrict const int* EToBM,
@restrict const int* EToB,
@restrict const int* EToBV,
@restrict const dfloat* W,
@restrict const dfloat* U,
@restrict dfloat* UH)
......@@ -277,6 +278,7 @@
const dlong sid = e * p_Nfaces * p_Nfp + n;
const dlong idM = vmapM[sk];
const dlong bcType = EToB[f + p_Nfaces * e];
const dlong bcTypeV = EToBV[idM];
if(bcType > 0) {
bc.id = EToBM[f + p_Nfaces * e];
......@@ -293,26 +295,30 @@
bc.u = U[idM + 0 * offset];
bc.v = U[idM + 1 * offset];
bc.w = U[idM + 2 * offset];
//bc.pP =
bc.wrk = W;
bc.fieldOffset = offset;
if(bcType == 1) {
bc.u = 0.f;
bc.v = 0.f;
bc.w = 0.f;
}else if(bcType == 2) {
UH[idM + 0 * offset] = 0.0;
UH[idM + 1 * offset] = 0.0;
UH[idM + 2 * offset] = 0.0;
}else if(bcType == 2 && bcTypeV == bcType) {
velocityDirichletConditions(&bc);
}else if(bcType == 4) {
bc.u = 0.f; // vP = vM ; wP = wM;
}else if(bcType == 5) {
bc.v = 0.f; // uP = vM ; wP = wM;
}else if(bcType == 6) {
bc.w = 0.f; // vP = vM ; uP = uM;
UH[idM + 0 * offset] = bc.u;
UH[idM + 1 * offset] = bc.v;
UH[idM + 2 * offset] = bc.w;
}else if(bcType == 3 && bcTypeV == bcType) {
UH[idM + 0 * offset] = bc.u;
UH[idM + 1 * offset] = bc.v;
UH[idM + 2 * offset] = bc.w;
}else if(bcType == 4 && bcTypeV == bcType) {
UH[idM + 0 * offset] = 0.0;
}else if(bcType == 5 && bcTypeV == bcType) {
UH[idM + 1 * offset] = 0.0;
}else if(bcType == 6 && bcTypeV == bcType) {
UH[idM + 2 * offset] = 0.0;
}
UH[idM + 0 * offset] = bc.u;
UH[idM + 1 * offset] = bc.v;
UH[idM + 2 * offset] = bc.w;
}
}
}
......
......@@ -57,8 +57,6 @@ occa::memory cdsSolve(const int is, cds_t* cds, dfloat time)
cds->o_mapB[is],
*(cds->o_usrwrk),
cds->o_wrk1);
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[is].compareArgs("SCALAR INITIAL GUESS DEFAULT", "EXTRAPOLATION")) {
cds->o_wrk0.copyFrom(cds->o_Se, cds->Ntotal * sizeof(dfloat), 0, is * cds->fieldOffset * sizeof(dfloat));
......
......@@ -41,7 +41,6 @@ void ellipticAx(elliptic_t* elliptic,
const int continuous = options.compareArgs("DISCRETIZATION", "CONTINUOUS");
const int serial = options.compareArgs("THREAD MODEL", "SERIAL");
const int ipdg = options.compareArgs("DISCRETIZATION", "IPDG");
const int mapType = (elliptic->elementType == HEXAHEDRA &&
options.compareArgs("ELEMENT MAP", "TRILINEAR")) ? 1:0;
const int integrationType = (elliptic->elementType == HEXAHEDRA &&
......
......@@ -33,33 +33,32 @@ int ellipticSolve(elliptic_t* elliptic,
mesh_t* mesh = elliptic->mesh;
setupAide options = elliptic->options;
int Niter = 0;
int maxIter = 1000;
dfloat tol = 1e-6;
options.getArgs("MAXIMUM ITERATIONS", maxIter);
dfloat tol = 1e-6;
options.getArgs("SOLVER TOLERANCE", tol);
elliptic->resNormFactor = 1 / (elliptic->Nfields * mesh->volume);
if(elliptic->var_coeff && options.compareArgs("PRECONDITIONER", "JACOBI"))
ellipticUpdateJacobi(elliptic);
if(options.compareArgs("RESIDUAL PROJECTION","TRUE"))
elliptic->o_x0.copyFrom(o_x, elliptic->Nfields * elliptic->Ntotal * sizeof(dfloat));
// compute initial residual
ellipticOperator(elliptic, o_x, elliptic->o_Ap, dfloatString);
// compute initial residual r = rhs - x0
ellipticAx(elliptic, mesh->NglobalGatherElements, mesh->o_globalGatherElementList, o_x, elliptic->o_Ap, dfloatString);
ellipticAx(elliptic, mesh->NlocalGatherElements, mesh->o_localGatherElementList, o_x, elliptic->o_Ap, dfloatString);
ellipticScaledAdd(elliptic, -1.f, elliptic->o_Ap, 1.f, o_r);
if (elliptic->Nmasked) mesh->maskKernel(elliptic->Nmasked, elliptic->o_maskIds, o_r);
if(elliptic->allNeumann)
ellipticZeroMean(elliptic, o_r);
if(elliptic->allNeumann) ellipticZeroMean(elliptic, o_r);
oogs::startFinish(o_r, elliptic->Nfields, elliptic->Ntotal, ogsDfloat, ogsAdd, elliptic->oogs);
if(elliptic->Nmasked) mesh->maskKernel(elliptic->Nmasked, elliptic->o_maskIds, o_r);
if(options.compareArgs("RESIDUAL PROJECTION","TRUE")) {
timer::tic("pre",1);
elliptic->o_x0.copyFrom(o_x, elliptic->Nfields * elliptic->Ntotal * sizeof(dfloat));
elliptic->residualProjection->pre(o_r);
timer::toc("pre");
}
dlong Niter;
if(!options.compareArgs("KRYLOV SOLVER", "NONBLOCKING")) {
Niter = pcg (elliptic, o_r, o_x, tol, maxIter);
}else{
......
......@@ -25,9 +25,8 @@
*/
#include "elliptic.h"
//#include "ogsInterface.h"
#define USE_WEIGHTED 1
#define USE_WEIGHTED 0
void ellipticZeroMean(elliptic_t* elliptic, occa::memory &o_q)
{
......
......@@ -252,10 +252,9 @@ void printRuntimeStatistics()
static void dryRun(setupAide &options, int npTarget)
{
if (rank == 0)
cout << "performing dry-run to jit-compile for >"
<< npTarget
<< " MPI tasks ...\n" << endl;
cout << "performing dry-run to jit-compile for >"
<< npTarget
<< " MPI tasks ...\n" << endl;
options.setArgs("NP TARGET", std::to_string(npTarget));
options.setArgs("BUILD ONLY", "TRUE");
......@@ -264,7 +263,7 @@ static void dryRun(setupAide &options, int npTarget)
string udfFile;
options.getArgs("UDF FILE", udfFile);
if (!udfFile.empty()) {
if(rank == 0) udfBuild(udfFile.c_str());
udfBuild(udfFile.c_str());
MPI_Barrier(comm);
*(void**)(&udf.loadKernels) = udfLoadFunction("UDF_LoadKernels",0);
*(void**)(&udf.setup0) = udfLoadFunction("UDF_Setup0",0);
......@@ -275,7 +274,7 @@ static void dryRun(setupAide &options, int npTarget)
// init solver
nrsSetup(comm, device, options, nrs);
if (rank == 0) cout << "\nBuild successful." << endl;
cout << "\nBuild successful." << endl;
}
static void setOUDF(setupAide &options)
......
......@@ -35,6 +35,7 @@ occa::memory pressureSolve(nrs_t* nrs, dfloat time)
mesh->o_vmapM,
mesh->o_EToB,
nrs->o_EToB,
nrs->o_VmapB,
nrs->o_usrwrk,
nrs->o_U,
nrs->o_wrk7);
......@@ -151,8 +152,6 @@ occa::memory pressureSolve(nrs_t* nrs, dfloat time)
nrs->o_U,
nrs->o_wrk3);
oogs::startFinish(nrs->o_wrk3, 1, 0, ogsDfloat, ogsAdd, nrs->gsh);
nrs->o_wrk1.copyFrom(nrs->o_P, nrs->Ntotal * sizeof(dfloat));
nrs->NiterP = ellipticSolve(nrs->pSolver, nrs->o_wrk3, nrs->o_wrk1);
......@@ -239,8 +238,6 @@ occa::memory velocitySolve(nrs_t* nrs, dfloat time)
nrs->o_rho,
nrs->o_wrk3);
oogs::startFinish(nrs->o_wrk3, nrs->NVfields, nrs->fieldOffset,ogsDfloat, ogsAdd, nrs->gsh);
if(nrs->options.compareArgs("VELOCITY INITIAL GUESS DEFAULT", "EXTRAPOLATION")) {
nrs->o_wrk0.copyFrom(nrs->o_Ue, nrs->NVfields * nrs->fieldOffset * sizeof(dfloat));
if (nrs->uvwSolver) {
......
......@@ -21,6 +21,8 @@ $(OBJDIR)/nekInterface.o: $(NEKRS_NEKINTERFACE_DIR)/nekInterface.f $(LIB)
libnekInterface: lib$(CASENAME).so
lib$(CASENAME).so: $(LIB) $(USRF) $(OBJDIR)/nekInterface.o
$(FC) $(FL2) -shared -o lib$(CASENAME).so $(USRF) $(OBJDIR)/nekInterface.o $(NI_LDFLAGS)
$(FC) $(FL2) -shared -o lib$(CASENAME).so $(USRF) $(OBJDIR)/nekInterface.o $(NI_LDFLAGS)
.NOTPARALLEL:
.PHONY: libnekInterface
......@@ -586,9 +586,6 @@ int nek_setup(MPI_Comm c, setupAide &options_in, nrs_t* nrs_in)
dfloat lambda;
options->getArgs("SCALAR00 DIFFUSIVITY", lambda);
if(!bcInPar)
(*nek_setup_ptr)(&nek_comm, (char*)cwd.c_str(), (char*)casename.c_str(),
&flow, &nscal, &nBcRead, &meshPartType,
&rho, &mue, &rhoCp, &lambda,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment