is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| Cal_SolutionErrorRatio.cpp (getdp-3.4.0-source.tgz) | : | Cal_SolutionErrorRatio.cpp (getdp-3.5.0-source.tgz) |
|---|
|
| // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege | | // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege |
| // | | // |
| // See the LICENSE.txt file for license information. Please report all | | // See the LICENSE.txt file for license information. Please report all |
| // issues on https://gitlab.onelab.info/getdp/getdp/issues. | | // issues on https://gitlab.onelab.info/getdp/getdp/issues. |
| // | | // |
| // Contributor(s): | | // Contributor(s): |
| // Michael Asam | | // Michael Asam |
| | |
| #include <stdio.h> | | #include <stdio.h> |
| #include <limits> | | #include <limits> |
| #include <math.h> | | #include <math.h> |
| #include "ProData.h" | | #include "ProData.h" |
| #include "DofData.h" | | #include "DofData.h" |
| #include "SolvingOperations.h" | | #include "SolvingOperations.h" |
| #include "Message.h" | | #include "Message.h" |
| | |
|
| void Cal_SolutionErrorRatio(gVector *dx, gVector *x, | | void Cal_SolutionErrorRatio(gVector *dx, gVector *x, double reltol, |
| double reltol, double abstol, | | double abstol, int NormType, double *ErrorRatio) |
| int NormType, double *ErrorRatio) | | |
| { | | { |
|
| int xLength; | | int xLength; |
| double AbsVal_x, AbsVal_dx, ImagVal_x, ImagVal_dx; | | double AbsVal_x, AbsVal_dx, ImagVal_x, ImagVal_dx; |
| double *ErrorRatioVec; | | double *ErrorRatioVec; |
| bool Is_NaN_or_Inf; | | bool Is_NaN_or_Inf; |
| | |
| LinAlg_GetVectorSize(dx, &xLength); | | LinAlg_GetVectorSize(dx, &xLength); |
| ErrorRatioVec = new double[xLength]; | | ErrorRatioVec = new double[xLength]; |
| | |
| *ErrorRatio = 0.; | | *ErrorRatio = 0.; |
| Is_NaN_or_Inf = false; | | Is_NaN_or_Inf = false; |
| | |
|
| for (int i = 0; i < xLength; i++) { | | for(int i = 0; i < xLength; i++) { |
| if (gSCALAR_SIZE == 1) | | if(gSCALAR_SIZE == 1) { |
| { | | LinAlg_GetAbsDoubleInVector(&AbsVal_x, x, i); |
| LinAlg_GetAbsDoubleInVector(&AbsVal_x, x, i) ; | | LinAlg_GetAbsDoubleInVector(&AbsVal_dx, dx, i); |
| LinAlg_GetAbsDoubleInVector(&AbsVal_dx, dx, i) ; | | |
| } | | } |
|
| if (gSCALAR_SIZE == 2) | | if(gSCALAR_SIZE == 2) { |
| { | | |
| LinAlg_GetComplexInVector(&AbsVal_x, &ImagVal_x, x, i, -1); | | LinAlg_GetComplexInVector(&AbsVal_x, &ImagVal_x, x, i, -1); |
| LinAlg_GetComplexInVector(&AbsVal_dx, &ImagVal_dx, dx, i, -1); | | LinAlg_GetComplexInVector(&AbsVal_dx, &ImagVal_dx, dx, i, -1); |
|
| AbsVal_x = sqrt( AbsVal_x*AbsVal_x + ImagVal_x*ImagVal_x); | | AbsVal_x = sqrt(AbsVal_x * AbsVal_x + ImagVal_x * ImagVal_x); |
| AbsVal_dx = sqrt( AbsVal_dx*AbsVal_dx + ImagVal_dx*ImagVal_dx); | | AbsVal_dx = sqrt(AbsVal_dx * AbsVal_dx + ImagVal_dx * ImagVal_dx); |
| } | | } |
| | |
| ErrorRatioVec[i] = AbsVal_dx / (abstol + reltol * AbsVal_x); | | ErrorRatioVec[i] = AbsVal_dx / (abstol + reltol * AbsVal_x); |
| | |
|
| if ( ErrorRatioVec[i] != ErrorRatioVec[i] || // So | | if(ErrorRatioVec[i] != ErrorRatioVec[i] || // Solution is NaN |
| lution is NaN | | ErrorRatioVec[i] == |
| ErrorRatioVec[i] == -std::numeric_limits<double>::infinity() || // So | | -std::numeric_limits<double>::infinity() || // Solution is -Inf |
| lution is -Inf | | ErrorRatioVec[i] == |
| ErrorRatioVec[i] == std::numeric_limits<double>::infinity() ) // So | | std::numeric_limits<double>::infinity()) // Solution is Inf |
| lution is Inf | | |
| Is_NaN_or_Inf = true; | | Is_NaN_or_Inf = true; |
| } | | } |
| | |
|
| if ( Is_NaN_or_Inf ) { | | if(Is_NaN_or_Inf) { |
| //Message::Error("No valid solution found (NaN or Inf)!"); // kj+++ | | |
| Message::Warning("No valid solution found (NaN or Inf)!"); | | Message::Warning("No valid solution found (NaN or Inf)!"); |
| *ErrorRatio = std::numeric_limits<double>::quiet_NaN(); | | *ErrorRatio = std::numeric_limits<double>::quiet_NaN(); |
| } | | } |
|
| else{ | | else { |
| // Calculating the norm of the error ratio vector | | // Calculating the norm of the error ratio vector |
|
| switch (NormType){ | | switch(NormType) { |
| case LINFNORM: | | case LINFNORM: |
|
| for (int i = 0; i < xLength; i++) { | | for(int i = 0; i < xLength; i++) { |
| if (ErrorRatioVec[i] > *ErrorRatio) | | if(ErrorRatioVec[i] > *ErrorRatio) *ErrorRatio = ErrorRatioVec[i]; |
| *ErrorRatio = ErrorRatioVec[i]; | | |
| } | | } |
| break; | | break; |
| | |
| case L1NORM: | | case L1NORM: |
|
| for (int i = 0; i < xLength; i++) { | | for(int i = 0; i < xLength; i++) { *ErrorRatio += ErrorRatioVec[i]; } |
| *ErrorRatio += ErrorRatioVec[i]; | | |
| } | | |
| break; | | break; |
| | |
| case MEANL1NORM: | | case MEANL1NORM: |
|
| for (int i = 0; i < xLength; i++) { | | for(int i = 0; i < xLength; i++) { *ErrorRatio += ErrorRatioVec[i]; } |
| *ErrorRatio += ErrorRatioVec[i]; | | |
| } | | |
| *ErrorRatio /= xLength; | | *ErrorRatio /= xLength; |
| break; | | break; |
| | |
| case L2NORM: | | case L2NORM: |
|
| for (int i = 0; i < xLength; i++) { | | for(int i = 0; i < xLength; i++) { |
| *ErrorRatio += ErrorRatioVec[i] * ErrorRatioVec[i]; | | *ErrorRatio += ErrorRatioVec[i] * ErrorRatioVec[i]; |
| } | | } |
| *ErrorRatio = sqrt(*ErrorRatio); | | *ErrorRatio = sqrt(*ErrorRatio); |
| break; | | break; |
| | |
| case MEANL2NORM: | | case MEANL2NORM: |
|
| for (int i = 0; i < xLength; i++) { | | for(int i = 0; i < xLength; i++) { |
| *ErrorRatio += ErrorRatioVec[i] * ErrorRatioVec[i]; | | *ErrorRatio += ErrorRatioVec[i] * ErrorRatioVec[i]; |
| } | | } |
| *ErrorRatio = sqrt(*ErrorRatio / xLength); | | *ErrorRatio = sqrt(*ErrorRatio / xLength); |
| break; | | break; |
| | |
| default: | | default: |
| Message::Error("Wrong error norm in Cal_SolutionErrorRatio"); | | Message::Error("Wrong error norm in Cal_SolutionErrorRatio"); |
| break; | | break; |
| } | | } |
| } | | } |
| | |
|
| delete [] ErrorRatioVec; | | delete[] ErrorRatioVec; |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l _ S o l u t i o n E r r o r */ | | /* C a l _ S o l u t i o n E r r o r */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
| void Cal_SolutionError(gVector *dx, gVector *x, int diff, double *MeanError) | | void Cal_SolutionError(gVector *dx, gVector *x, int diff, double *MeanError) |
| { | | { |
| // This is not a very good implementation: it should be replaced with | | // This is not a very good implementation: it should be replaced with |
| // Cal_SolutionErrorRatio above | | // Cal_SolutionErrorRatio above |
| | |
|
| int i, n; | | int i, n; |
| double valx, valdx, valxi = 0., valdxi = 0.,errsqr = 0., xmoy = 0., dxmoy = 0. | | double valx, valdx, valxi = 0., valdxi = 0., errsqr = 0., xmoy = 0., |
| ; | | dxmoy = 0.; |
| double tol, nvalx, nvaldx ; | | double tol, nvalx, nvaldx; |
| //bool Is_NaN_or_Inf=0; // kj+++ | | |
| | |
| LinAlg_GetVectorSize(dx, &n); | | LinAlg_GetVectorSize(dx, &n); |
| | |
|
| if (gSCALAR_SIZE == 1) | | if(gSCALAR_SIZE == 1) |
| for (i=0 ; i<n ; i++) { | | for(i = 0; i < n; i++) { |
| LinAlg_GetAbsDoubleInVector(&valx, x, i) ; | | LinAlg_GetAbsDoubleInVector(&valx, x, i); |
| LinAlg_GetAbsDoubleInVector(&valdx, dx, i) ; | | LinAlg_GetAbsDoubleInVector(&valdx, dx, i); |
| xmoy += valx ; | | xmoy += valx; |
| if(diff) dxmoy += (valdx-valx) ; | | if(diff) |
| else dxmoy += valdx ; | | dxmoy += (valdx - valx); |
| | else |
| | dxmoy += valdx; |
| } | | } |
|
| if (gSCALAR_SIZE == 2) | | if(gSCALAR_SIZE == 2) |
| for (i=0 ; i<n ; i++) { | | for(i = 0; i < n; i++) { |
| LinAlg_GetComplexInVector(&valx, &valxi, x, i, -1); | | LinAlg_GetComplexInVector(&valx, &valxi, x, i, -1); |
| LinAlg_GetComplexInVector(&valdx, &valdxi, dx, i, -1); | | LinAlg_GetComplexInVector(&valdx, &valdxi, dx, i, -1); |
|
| xmoy += sqrt(valx*valx+valxi*valxi) ; | | xmoy += sqrt(valx * valx + valxi * valxi); |
| if(diff) dxmoy += sqrt((valdx-valx)*(valdx-valx)+(valdxi-valxi)*(valdxi-va | | if(diff) |
| lxi)) ; | | dxmoy += sqrt((valdx - valx) * (valdx - valx) + |
| else dxmoy += sqrt(valdx*valdx + valdxi*valdxi) ; | | (valdxi - valxi) * (valdxi - valxi)); |
| | else |
| | dxmoy += sqrt(valdx * valdx + valdxi * valdxi); |
| } | | } |
| | |
|
| xmoy /= (double)n ; | | xmoy /= (double)n; |
| dxmoy /= (double)n ; | | dxmoy /= (double)n; |
| | |
|
| if (xmoy > 1.e-30) { | | if(xmoy > 1.e-30) { |
| tol = xmoy*1.e-10 ; | | tol = xmoy * 1.e-10; |
| //tol = 1e200; // kj+++ | | if(gSCALAR_SIZE == 1) |
| //tol=(tol<1e-3)?1e-3:tol; // kj+++ | | for(i = 0; i < n; i++) { |
| if (gSCALAR_SIZE == 1) | | LinAlg_GetAbsDoubleInVector(&valx, x, i); |
| for (i=0 ; i<n ; i++){ | | LinAlg_GetAbsDoubleInVector(&valdx, dx, i); |
| LinAlg_GetAbsDoubleInVector(&valx, x, i) ; | | if(diff) { |
| LinAlg_GetAbsDoubleInVector(&valdx, dx, i) ; | | if(valx > tol) |
| if(diff){ | | errsqr += fabs(valdx - valx) / valx; |
| if (valx > tol) errsqr += fabs(valdx-valx)/valx ; | | else |
| else errsqr += fabs(valdx-valx) ; | | errsqr += fabs(valdx - valx); |
| } | | } |
|
| else{ | | else { |
| if (valx > tol) errsqr += valdx/valx ; | | if(valx > tol) |
| else errsqr += valdx ; | | errsqr += valdx / valx; |
| | else |
| | errsqr += valdx; |
| } | | } |
| } | | } |
| | |
|
| if (gSCALAR_SIZE == 2) | | if(gSCALAR_SIZE == 2) |
| for (i=0 ; i<n ; i++) { | | for(i = 0; i < n; i++) { |
| LinAlg_GetComplexInVector(&valx, &valxi, x, i, -1); | | LinAlg_GetComplexInVector(&valx, &valxi, x, i, -1); |
| LinAlg_GetComplexInVector(&valdx, &valdxi, dx, i, -1); | | LinAlg_GetComplexInVector(&valdx, &valdxi, dx, i, -1); |
|
| nvalx = sqrt(valx*valx+valxi*valxi) ; | | nvalx = sqrt(valx * valx + valxi * valxi); |
| nvaldx = sqrt(valdx*valdx+valdxi*valdxi) ; | | nvaldx = sqrt(valdx * valdx + valdxi * valdxi); |
| if(diff){ | | if(diff) { |
| if (nvalx > tol) | | if(nvalx > tol) |
| errsqr += sqrt((valdx-valx)*(valdx-valx)+(valdxi-valxi)*(valdxi-valx | | errsqr += sqrt((valdx - valx) * (valdx - valx) + |
| i))/nvalx ; | | (valdxi - valxi) * (valdxi - valxi)) / |
| | nvalx; |
| else | | else |
|
| errsqr += sqrt((valdx-valx)*(valdx-valx)+(valdxi-valxi)*(valdxi-valx | | errsqr += sqrt((valdx - valx) * (valdx - valx) + |
| i)); | | (valdxi - valxi) * (valdxi - valxi)); |
| } | | } |
|
| else{ | | else { |
| if (nvalx > tol) errsqr += nvaldx/nvalx ; | | if(nvalx > tol) |
| else errsqr += nvaldx ; | | errsqr += nvaldx / nvalx; |
| | else |
| | errsqr += nvaldx; |
| } | | } |
| } | | } |
| | |
|
| *MeanError = errsqr/(double)n ; | | *MeanError = errsqr / (double)n; |
| } | | } |
|
| else{ | | else { |
| if (dxmoy > 1.e-30) | | if(dxmoy > 1.e-30) |
| *MeanError = 1. ; | | *MeanError = 1.; |
| else | | else |
|
| *MeanError = 0. ; | | *MeanError = 0.; |
| } | | |
| | |
| /* // kj+++ | | |
| if ( *MeanError != *MeanError || // Solution i | | |
| s NaN | | |
| *MeanError == -std::numeric_limits<double>::infinity() || // Solution i | | |
| s -Inf | | |
| *MeanError == std::numeric_limits<double>::infinity() ) // Solution i | | |
| s Inf | | |
| Is_NaN_or_Inf = true; | | |
| | |
| if ( Is_NaN_or_Inf ) { | | |
| Message::Error("No valid solution found (NaN or Inf)!"); | | |
| *MeanError = std::numeric_limits<double>::quiet_NaN(); | | |
| } | | } |
|
| */ | | |
| } | | } |
| | | |
| End of changes. 31 change blocks. |
|---|
| 106 lines changed or deleted | | 89 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 