Commit 29a99c87 authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard
Browse files

KernelGen Nonlinear STP LoG - test advanced stop criterion (not very good if zero values expected)

parent da5a29f8
......@@ -102,6 +102,7 @@ class Controller:
"useVectPDE" : args["useVectPDE"],
"useAoSoA2" : args["useAoSoA2"],
"predictorRecompute" : args["predictorRecompute"],
"advancedStopCriterion" : False, #TODO JMG put as proper toolkit arg
"initialGuess" : "mixedPicard" #TODO JMG put as proper toolkit arg
#"initialGuess" : "default" #TODO JMG put as proper toolkit arg
})
......
......@@ -189,7 +189,9 @@ int {{codeNamespace}}::fusedSpaceTimePredictorVolumeIntegral(
{% endif %}{# useCERKGuess #}
// 2. Discrete Picard iterations
constexpr int MaxIterations = {% if useCERKGuess %}{% if nDof-3 <= 1 %}1; //cannot be lower than 1{% else %}{{nDof-3}}; //nDof-3{% endif %}{% else %}{{2*nDof+1}};{% endif %}
//constexpr int MaxIterations = {% if useCERKGuess %}{% if nDof-3 <= 1 %}1; //cannot be lower than 1{% else %}{{nDof-3}}; //nDof-3{% endif %}{% else %}{{2*nDof+1}};{% endif %}
constexpr int MaxIterations = {{2*nDof+1}};
int iter = 0;
for (; iter < MaxIterations; iter++) {
......@@ -308,23 +310,45 @@ int {{codeNamespace}}::fusedSpaceTimePredictorVolumeIntegral(
// 3. Multiply with (K1)^(-1) to get the discrete time integral of the
// discrete Picard iteration
double sq_res = 0.0;
{% if advancedStopCriterion %}
double max_approx = 0.;
// for (int xyzt = 0; xyzt < {{(nDof**nDim)*nDof}}; xyzt++) {
// for(int n=0; n<{{nVar}}; n++) {
// max_approx += lQi[n+{{nDataPad}}*xyzt]*lQi[n+{{nDataPad}}*xyzt];
// }
// }
// max_approx = sqrt(max_approx)/{{(nDof**nDim)*nDof*nVar}};
max_approx = std::fabs(*std::max_element(lQi, lQi+{{(nDof**nDim)*nDof*nDataPad}}, [](float a, float b){return (std::fabs(a) < std::fabs(b));}));
double machine_precision = max_approx*1e-10; // used to avoid rounding error relative error to prevent convergence, assume all variables are correlated (and thus value below avg variable * 10^-10 is a rounding error)
{% endif %}
for (int xyz = 0; xyz < {{nDof**nDim}}; xyz++) {
{{ m.matmul_legacy('lqi', 'rhs', 's_m_QSlice', 'new_lQi_slice',nVarPad~'*xyz', '0', '0', trueB='iK1_T', trueAlpha='iweights3[xyz]') | indent(6) }}{##}
for(int t = 0; t < {{nDof}}; t++) {
for(int n=0; n<{{nVar}}; n++) { //only copy and change the variables, skip parameters
{% if advancedStopCriterion %}
sq_res += ((new_lQi_slice[n+{{nVarPad}}*t] - lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)]) * (new_lQi_slice[n+{{nVarPad}}*t] - lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)]))/((lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)]*lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)])+machine_precision);
{% else %}
sq_res += (new_lQi_slice[n+{{nVarPad}}*t] - lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)]) * (new_lQi_slice[n+{{nVarPad}}*t] - lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)]);
{% endif %}
lQi[n+{{nDataPad}}*(xyz+{{nDof**nDim}}*t)] = new_lQi_slice[n+{{nVarPad}}*t];
}
}
}
// 4. Exit condition
{% if advancedStopCriterion %}
constexpr double tol2 = 1e-11 * 1e-11 * {{(nDof**nDim)*nDof*nVar}};
{% else %}
constexpr double tol2 = 1e-7 * 1e-7;
{% endif %}
if (sq_res < tol2) {
break;
}
} // end iter
//if (iter == MaxIterations){printf("MaxIter reached \n");}
//***********************
//****** Predictor ******
//***********************
......
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