is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| Operation_TimeLoopAdaptive.cpp (getdp-3.4.0-source.tgz) | : | Operation_TimeLoopAdaptive.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 <stdlib.h> | | #include <stdlib.h> |
| #include <string.h> | | #include <string.h> |
| | |
| skipping to change at line 28 | | skipping to change at line 28 |
|---|
| #include "Message.h" | | #include "Message.h" |
| #include "MallocUtils.h" | | #include "MallocUtils.h" |
| #include "Legendre.h" | | #include "Legendre.h" |
| | |
| #if !defined(F77NAME) | | #if !defined(F77NAME) |
| #define F77NAME(x) (x##_) | | #define F77NAME(x) (x##_) |
| #endif | | #endif |
| | |
| #if defined(HAVE_LAPACK) | | #if defined(HAVE_LAPACK) |
| extern "C" { | | extern "C" { |
|
| void F77NAME(dgesv)(int *N, int *nrhs, double *A, int *lda, int *ipiv, | | void F77NAME(dgesv)(int *N, int *nrhs, double *A, int *lda, int *ipiv, |
| double *b, int *ldb, int *info); | | double *b, int *ldb, int *info); |
| } | | } |
| #endif | | #endif |
| | |
| extern struct CurrentData Current; | | extern struct CurrentData Current; |
| extern int Flag_IterativeLoopConverged; | | extern int Flag_IterativeLoopConverged; |
| extern int Flag_RESTART; | | extern int Flag_RESTART; |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l c I n t e g r a t i o n C o e f f i c i e n t s */ | | /* C a l c I n t e g r a t i o n C o e f f i c i e n t s */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
| #if !defined(HAVE_LAPACK) | | #if !defined(HAVE_LAPACK) |
| | |
|
| void CalcIntegrationCoefficients(Resolution *Resolution_P, | | void CalcIntegrationCoefficients(Resolution *Resolution_P, DofData *DofData_P0, |
| DofData *DofData_P0, | | List_T *TLAsystems_L, List_T *TLAPostOp_L, |
| List_T *TLAsystems_L, | | int Order) |
| List_T *TLAPostOp_L, | | |
| int Order) | | |
| { | | { |
| Message::Error("TimeLoopAdaptive requires Lapack"); | | Message::Error("TimeLoopAdaptive requires Lapack"); |
| } | | } |
| | |
| #else | | #else |
| | |
|
| void CalcIntegrationCoefficients(Resolution *Resolution_P, | | void CalcIntegrationCoefficients(Resolution *Resolution_P, DofData *DofData_P0, |
| DofData *DofData_P0, | | List_T *TLAsystems_L, List_T *TLAPostOp_L, |
| List_T *TLAsystems_L, | | int Order) |
| List_T *TLAPostOp_L, | | |
| int Order) | | |
| { | | { |
|
| DefineSystem *System_P=NULL; | | DefineSystem *System_P = NULL; |
| DofData *DofData_P=NULL; | | DofData *DofData_P = NULL; |
| Solution *Solution_P; | | Solution *Solution_P; |
| TimeLoopAdaptiveSystem TLAsystem; | | TimeLoopAdaptiveSystem TLAsystem; |
|
| List_T *Solutions_L=NULL; | | List_T *Solutions_L = NULL; |
| PostOpSolutions *PostOpSolutions_P0; | | PostOpSolutions *PostOpSolutions_P0; |
| int j, NbrOfRows, NbrOfSolutions=0; | | int j, NbrOfRows, NbrOfSolutions = 0; |
| int Info, Pivots[7], NbrOfRightHandSides=1; | | int Info, Pivots[7], NbrOfRightHandSides = 1; |
| double t[8], temp; | | double t[8], temp; |
| double A[49], b[7]; | | double A[49], b[7]; |
| bool RecomputeTimeStep; | | bool RecomputeTimeStep; |
| | |
| // Initialization | | // Initialization |
|
| for (int i=0; i<7; i++) | | for(int i = 0; i < 7; i++) Current.aPredCoeff[i] = 0.0; |
| Current.aPredCoeff[i] = 0.0; | | for(int i = 0; i < 6; i++) Current.aCorrCoeff[i] = 0.0; |
| for (int i=0; i<6; i++) | | |
| Current.aCorrCoeff[i] = 0.0; | | |
| Current.bCorrCoeff = 0.0; | | Current.bCorrCoeff = 0.0; |
| Current.PredErrorConst = 0.0; | | Current.PredErrorConst = 0.0; |
| Current.CorrErrorConst = 0.0; | | Current.CorrErrorConst = 0.0; |
| | |
|
| if (Order < 1 || Order > 6) | | if(Order < 1 || Order > 6) |
| Message::Error("Order has to be in the range 1 .. 6"); | | Message::Error("Order has to be in the range 1 .. 6"); |
| | |
| // First get the past time points | | // First get the past time points |
| // ------------------------------ | | // ------------------------------ |
| | |
|
| if (List_Nbr(TLAsystems_L)==0 && List_Nbr(TLAPostOp_L)==0) { | | if(List_Nbr(TLAsystems_L) == 0 && List_Nbr(TLAPostOp_L) == 0) { |
| Message::Error("Neither systems nor PostOperations are specified " | | Message::Error("Neither systems nor PostOperations are specified " |
|
| "for TimeLoopAdaptive"); | | "for TimeLoopAdaptive"); |
| } | | } |
|
| if (List_Nbr(TLAsystems_L)) { | | if(List_Nbr(TLAsystems_L)) { |
| List_Read(TLAsystems_L, 0, &TLAsystem); | | List_Read(TLAsystems_L, 0, &TLAsystem); |
|
| System_P = (DefineSystem*)List_Pointer(Resolution_P->DefineSystem, | | System_P = (DefineSystem *)List_Pointer(Resolution_P->DefineSystem, |
| TLAsystem.SystemIndex); | | TLAsystem.SystemIndex); |
| DofData_P = DofData_P0 + TLAsystem.SystemIndex; | | DofData_P = DofData_P0 + TLAsystem.SystemIndex; |
| Solutions_L = DofData_P->Solutions; | | Solutions_L = DofData_P->Solutions; |
| NbrOfSolutions = List_Nbr(Solutions_L); | | NbrOfSolutions = List_Nbr(Solutions_L); |
|
| if (!NbrOfSolutions) | | if(!NbrOfSolutions) |
| Message::Error("No initial solution for system %s", System_P->Name); | | Message::Error("No initial solution for system %s", System_P->Name); |
|
| if (NbrOfSolutions <= Order && Order > 1) | | if(NbrOfSolutions <= Order && Order > 1) |
| Message::Error("Too few past solutions for system %s", System_P->Name); | | Message::Error("Too few past solutions for system %s", System_P->Name); |
| } | | } |
|
| if (List_Nbr(TLAPostOp_L)) { | | if(List_Nbr(TLAPostOp_L)) { |
| PostOpSolutions_P0 = (PostOpSolutions *)List_Pointer(Current.PostOpData_L, 0 | | PostOpSolutions_P0 = |
| ); | | (PostOpSolutions *)List_Pointer(Current.PostOpData_L, 0); |
| Solutions_L = PostOpSolutions_P0->Solutions_L; | | Solutions_L = PostOpSolutions_P0->Solutions_L; |
| NbrOfSolutions = List_Nbr(Solutions_L); | | NbrOfSolutions = List_Nbr(Solutions_L); |
|
| if (!NbrOfSolutions) | | if(!NbrOfSolutions) Message::Error("No initial PostOperations"); |
| Message::Error("No initial PostOperations"); | | if(NbrOfSolutions <= Order && Order > 1) |
| if (NbrOfSolutions <= Order && Order > 1) | | Message::Error("Too few past PostOperations results"); |
| Message::Error("Too few past PostOperations results"); | | |
| } | | } |
| | |
| // Set the predictor's and corrector's order | | // Set the predictor's and corrector's order |
| // ----------------------------------------- | | // ----------------------------------------- |
| | |
|
| Solution_P = (struct Solution*)List_Pointer(Solutions_L, NbrOfSolutions-1); | | Solution_P = (struct Solution *)List_Pointer(Solutions_L, NbrOfSolutions - 1); |
| | |
| // Check if we recompute actual TimeStep | | // Check if we recompute actual TimeStep |
| RecomputeTimeStep = (Solution_P->TimeStep == (int)Current.TimeStep); | | RecomputeTimeStep = (Solution_P->TimeStep == (int)Current.TimeStep); |
| | |
|
| if (NbrOfSolutions < (2 + (RecomputeTimeStep ? 1 : 0))){ | | if(NbrOfSolutions < (2 + (RecomputeTimeStep ? 1 : 0))) { |
| Current.PredOrder = 0; // For 1st TimeStep just copy the initial solution | | Current.PredOrder = 0; // For 1st TimeStep just copy the initial solution |
| Current.CorrOrder = 1; | | Current.CorrOrder = 1; |
| } | | } |
|
| else{ | | else { |
| Current.PredOrder = Order; | | Current.PredOrder = Order; |
| Current.CorrOrder = Order; | | Current.CorrOrder = Order; |
| } | | } |
| | |
| // Time values | | // Time values |
| // t_n+1 -> t[0] | | // t_n+1 -> t[0] |
| // t_n -> t[1] | | // t_n -> t[1] |
| // t_n-1 -> t[2] | | // t_n-1 -> t[2] |
| // ... | | // ... |
| // t_n-k -> t[k+1] k=Order | | // t_n-k -> t[k+1] k=Order |
| | |
| t[0] = Current.Time; | | t[0] = Current.Time; |
|
| for(int i=1; i <= Current.PredOrder+1; i++) { | | for(int i = 1; i <= Current.PredOrder + 1; i++) { |
| j = RecomputeTimeStep ? i+1 : i; | | j = RecomputeTimeStep ? i + 1 : i; |
| Solution_P = (struct Solution*)List_Pointer(Solutions_L, NbrOfSolutions-j); | | Solution_P = |
| | (struct Solution *)List_Pointer(Solutions_L, NbrOfSolutions - j); |
| t[i] = Solution_P->Time; | | t[i] = Solution_P->Time; |
| } | | } |
| | |
| // Calculation of predictor integration constants | | // Calculation of predictor integration constants |
| // ---------------------------------------------- | | // ---------------------------------------------- |
| | |
| /* The new solution is predicted by extrapolating the past solutions | | /* The new solution is predicted by extrapolating the past solutions |
| * by a polynom of order "PredOrder". The polynom coefficients | | * by a polynom of order "PredOrder". The polynom coefficients |
| * are calculated by solving a matrix equation A*coeff=b for the | | * are calculated by solving a matrix equation A*coeff=b for the |
| * exactness constraints. | | * exactness constraints. |
| | |
| skipping to change at line 162 | | skipping to change at line 157 |
|---|
| * | | * |
| * _ _ _ _ _ _ | | * _ _ _ _ _ _ |
| * | 1 1 1 1 | | a_0 | | 1 | | | * | 1 1 1 1 | | a_0 | | 1 | |
| * | (t_n)^1 (t_n-1)^1 (t_n-2)^1 (t_n-3)^1 | | a_1 | | (t_n+1)^1 | | | * | (t_n)^1 (t_n-1)^1 (t_n-2)^1 (t_n-3)^1 | | a_1 | | (t_n+1)^1 | |
| * | (t_n)^2 (t_n-1)^2 (t_n-2)^2 (t_n-3)^2 | * | a_2 | = | (t_n+1)^2 | | | * | (t_n)^2 (t_n-1)^2 (t_n-2)^2 (t_n-3)^2 | * | a_2 | = | (t_n+1)^2 | |
| * | (t_n)^3 (t_n-1)^3 (t_n-2)^3 (t_n-3)^3 | | a_3 | | (t_n+1)^3 | | | * | (t_n)^3 (t_n-1)^3 (t_n-2)^3 (t_n-3)^3 | | a_3 | | (t_n+1)^3 | |
| * |_ _| |_ _| |_ _| | | * |_ _| |_ _| |_ _| |
| * | | * |
| */ | | */ |
| | |
|
| if (Current.PredOrder == 0) | | if(Current.PredOrder == 0) |
| Current.aPredCoeff[0] = 1.0; | | Current.aPredCoeff[0] = 1.0; |
| else { | | else { |
| NbrOfRows = Current.PredOrder + 1; | | NbrOfRows = Current.PredOrder + 1; |
|
| for (int c=0; c <= Current.PredOrder; c++) { | | for(int c = 0; c <= Current.PredOrder; c++) { |
| A[0 + c*NbrOfRows] = 1.0; | | A[0 + c * NbrOfRows] = 1.0; |
| for (int r=1; r <= Current.PredOrder; r++) | | for(int r = 1; r <= Current.PredOrder; r++) |
| A[r + c*NbrOfRows] = pow(t[c+1],r); | | A[r + c * NbrOfRows] = pow(t[c + 1], r); |
| } | | } |
| | |
| b[0] = 1.0; | | b[0] = 1.0; |
|
| for (int r=1; r <= Current.PredOrder; r++) | | for(int r = 1; r <= Current.PredOrder; r++) b[r] = pow(t[0], r); |
| b[r] = pow(t[0],r); | | |
| | |
|
| F77NAME(dgesv)(&NbrOfRows, &NbrOfRightHandSides, A, &NbrOfRows, Pivots, b, & | | F77NAME(dgesv) |
| NbrOfRows, &Info); | | (&NbrOfRows, &NbrOfRightHandSides, A, &NbrOfRows, Pivots, b, &NbrOfRows, |
| if (Info != 0) | | &Info); |
| Message::Error("Can't calculate predictor coefficients for TimeLoopAdaptiv | | if(Info != 0) |
| e"); | | Message::Error( |
| | "Can't calculate predictor coefficients for TimeLoopAdaptive"); |
| | |
|
| for (int i=0; i<=Current.PredOrder; i++) | | for(int i = 0; i <= Current.PredOrder; i++) Current.aPredCoeff[i] = b[i]; |
| Current.aPredCoeff[i] = b[i]; | | |
| } | | } |
| | |
| // Calculation of corrector integration constants | | // Calculation of corrector integration constants |
| // ---------------------------------------------- | | // ---------------------------------------------- |
| | |
| /* | | /* |
| * The coefficients for the Gear method (BDF) are also | | * The coefficients for the Gear method (BDF) are also |
| * calculated by solving a matrix equation A*coeff=b for the | | * calculated by solving a matrix equation A*coeff=b for the |
| * exactness constraints. | | * exactness constraints. |
| * E.g. for CorrOder=3 we have: | | * E.g. for CorrOder=3 we have: |
| * | | * |
|
| * _ _ _ _ | | * _ _ _ _ _ _ |
| _ _ | | * | 1 1 1 0 | | a_0 | | |
| * | 1 1 1 0 | | a_0 | | | | * 1 | | (t_n)^1 (t_n-1)^1 (t_n-2)^1 1*(t_n+1 - t_n) | | |
| 1 | | | * a_1 | | (t_n+1)^1 | | (t_n)^2 (t_n-1)^2 (t_n-2)^2 2*(t_n+1 - |
| * | (t_n)^1 (t_n-1)^1 (t_n-2)^1 1*(t_n+1 - t_n) | | a_1 | | | | * t_n)*(t_n+1)^1 | * | a_2 | = | (t_n+1)^2 | | (t_n)^3 (t_n-1)^3 (t_n-2)^3 |
| (t_n+1)^1 | | | * 3*(t_n+1 - t_n)*(t_n+1)^2 | | b_-1 | | (t_n+1)^3 | |
| * | (t_n)^2 (t_n-1)^2 (t_n-2)^2 2*(t_n+1 - t_n)*(t_n+1)^1 | * | a_2 | = | | | * |_ _| |_ _| |_ |
| (t_n+1)^2 | | | * _| |
| * | (t_n)^3 (t_n-1)^3 (t_n-2)^3 3*(t_n+1 - t_n)*(t_n+1)^2 | | b_-1 | | | | |
| (t_n+1)^3 | | | |
| * |_ _| |_ _| | | | |
| _ _| | | |
| * | | * |
| */ | | */ |
| | |
|
| if (Current.TypeTime == TIME_GEAR) { | | if(Current.TypeTime == TIME_GEAR) { |
| NbrOfRows = Current.CorrOrder + 1; | | NbrOfRows = Current.CorrOrder + 1; |
|
| for (int c=0; c < Current.CorrOrder; c++) { | | for(int c = 0; c < Current.CorrOrder; c++) { |
| A[0 + c*NbrOfRows] = 1.0; | | A[0 + c * NbrOfRows] = 1.0; |
| for (int r=1; r <= Current.CorrOrder; r++) | | for(int r = 1; r <= Current.CorrOrder; r++) |
| A[r + c*NbrOfRows] = pow(t[c+1],r); | | A[r + c * NbrOfRows] = pow(t[c + 1], r); |
| } | | } |
| A[0 + (int)Current.CorrOrder*NbrOfRows] = 0.0; | | A[0 + (int)Current.CorrOrder * NbrOfRows] = 0.0; |
| A[1 + (int)Current.CorrOrder*NbrOfRows] = t[0]-t[1]; | | A[1 + (int)Current.CorrOrder * NbrOfRows] = t[0] - t[1]; |
| for (int r=2; r <= Current.CorrOrder; r++) | | for(int r = 2; r <= Current.CorrOrder; r++) |
| A[r + (int)Current.CorrOrder*NbrOfRows] = r*pow(t[0],r-1)*(t[0]-t[1]); | | A[r + (int)Current.CorrOrder * NbrOfRows] = |
| | r * pow(t[0], r - 1) * (t[0] - t[1]); |
| | |
| b[0] = 1.0; | | b[0] = 1.0; |
|
| for (int r=1; r <= Current.CorrOrder; r++) | | for(int r = 1; r <= Current.CorrOrder; r++) b[r] = pow(t[0], r); |
| b[r] = pow(t[0],r); | | |
| | |
|
| F77NAME(dgesv)(&NbrOfRows, &NbrOfRightHandSides, A, &NbrOfRows, Pivots, b, & | | F77NAME(dgesv) |
| NbrOfRows, &Info); | | (&NbrOfRows, &NbrOfRightHandSides, A, &NbrOfRows, Pivots, b, &NbrOfRows, |
| if (Info != 0) | | &Info); |
| Message::Error("Can't calculate corrector coefficients for TimeLoopAdaptiv | | if(Info != 0) |
| e"); | | Message::Error( |
| | "Can't calculate corrector coefficients for TimeLoopAdaptive"); |
| for (int i=0; i<Current.CorrOrder; i++) | | |
| Current.aCorrCoeff[i] = b[i]; | | for(int i = 0; i < Current.CorrOrder; i++) Current.aCorrCoeff[i] = b[i]; |
| Current.bCorrCoeff = b[(int)Current.CorrOrder]; | | Current.bCorrCoeff = b[(int)Current.CorrOrder]; |
| } | | } |
| | |
| // Calculation of predictor error constant | | // Calculation of predictor error constant |
| // ---------------------------------------------- | | // ---------------------------------------------- |
|
| for (int i=1; i<=Current.PredOrder; i++) | | for(int i = 1; i <= Current.PredOrder; i++) |
| Current.PredErrorConst += Current.aPredCoeff[i] * | | Current.PredErrorConst += |
| pow(t[1]-t[1+i], Current.PredOrder+1); | | Current.aPredCoeff[i] * pow(t[1] - t[1 + i], Current.PredOrder + 1); |
| Current.PredErrorConst *= pow(-1, Current.PredOrder) / | | Current.PredErrorConst *= |
| pow(t[0]-t[1], Current.PredOrder+1); | | pow(-1, Current.PredOrder) / pow(t[0] - t[1], Current.PredOrder + 1); |
| Current.PredErrorConst += 1; | | Current.PredErrorConst += 1; |
| Current.PredErrorConst /= Factorial(Current.PredOrder + 1); | | Current.PredErrorConst /= Factorial(Current.PredOrder + 1); |
| | |
| // Calculation of corrector error constant | | // Calculation of corrector error constant |
| // ---------------------------------------------- | | // ---------------------------------------------- |
|
| switch (Current.TypeTime) { | | switch(Current.TypeTime) { |
| case TIME_THETA: | | case TIME_THETA: |
| if (Current.CorrOrder == 1) | | if(Current.CorrOrder == 1) |
| Current.CorrErrorConst = -0.5; | | Current.CorrErrorConst = -0.5; |
| else if (Current.CorrOrder == 2) | | else if(Current.CorrOrder == 2) |
| Current.CorrErrorConst = -1./12.; | | Current.CorrErrorConst = -1. / 12.; |
| else | | else |
| Message::Error("Order %d not allowed for Theta scheme.", | | Message::Error("Order %d not allowed for Theta scheme.", |
| Current.CorrOrder); | | Current.CorrOrder); |
| break; | | break; |
| | | |
| case TIME_GEAR: | | case TIME_GEAR: |
| Current.CorrErrorConst = 1 / Factorial(Current.CorrOrder + 1); | | Current.CorrErrorConst = 1 / Factorial(Current.CorrOrder + 1); |
| | | |
| temp = 0.0; | | temp = 0.0; |
| for (int i=1; i<Current.CorrOrder; i++) | | for(int i = 1; i < Current.CorrOrder; i++) |
| temp += Current.aCorrCoeff[i] * pow(t[1]-t[1+i], Current.CorrOrder+1); | | temp += |
| temp *= pow(-1, Current.CorrOrder) / pow(t[0]-t[1], Current.CorrOrder+1); | | Current.aCorrCoeff[i] * pow(t[1] - t[1 + i], Current.CorrOrder + 1); |
| temp /= Factorial(Current.CorrOrder + 1); | | temp *= |
| Current.CorrErrorConst += temp; | | pow(-1, Current.CorrOrder) / pow(t[0] - t[1], Current.CorrOrder + 1); |
| | | temp /= Factorial(Current.CorrOrder + 1); |
| Current.CorrErrorConst -= Current.bCorrCoeff / | | Current.CorrErrorConst += temp; |
| Factorial(Current.CorrOrder); | | |
| break; | | Current.CorrErrorConst -= Current.bCorrCoeff / Factorial(Current.CorrOrder); |
| | | break; |
| default: | | |
| Message::Error("Unknown integration scheme for TimeLoopAdaptive"); | | default: |
| break; | | Message::Error("Unknown integration scheme for TimeLoopAdaptive"); |
| | break; |
| } | | } |
| } | | } |
| #endif | | #endif |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* P r e d i c t o r */ | | /* P r e d i c t o r */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Predictor(Resolution *Resolution_P, | | void Predictor(Resolution *Resolution_P, DofData *DofData_P0, |
| DofData *DofData_P0, | | List_T *TLAsystems_L, List_T *TLAPostOp_L, int Order, |
| List_T *TLAsystems_L, | | List_T *xPredicted_L, List_T *PostOpSolPredicted_L) |
| List_T *TLAPostOp_L, | | |
| int Order, | | |
| List_T *xPredicted_L, | | |
| List_T *PostOpSolPredicted_L) | | |
| { | | { |
|
| DefineSystem *System_P; | | DefineSystem *System_P; |
| DofData *DofData_P; | | DofData *DofData_P; |
| PostOpSolutions *PostOpSolutions_P; | | PostOpSolutions *PostOpSolutions_P; |
| Solution *Solution_P, *PastSolution_P, Solution_S; | | Solution *Solution_P, *PastSolution_P, Solution_S; |
| TimeLoopAdaptiveSystem TLAsystem; | | TimeLoopAdaptiveSystem TLAsystem; |
|
| gVector *xPredicted_P; | | gVector *xPredicted_P; |
| gVector *x_NminusJ_P; // past solution vector x_N- | | gVector *x_NminusJ_P; // past solution vector x_N-i |
| i | | gVector *PostOpSolPredicted_P; |
| gVector *PostOpSolPredicted_P; | | int TimeStep, NbrSolutions, PostOpSolLength; |
| int TimeStep, NbrSolutions, PostOpSolLength; | | |
| | |
| // Loop through all given systems | | // Loop through all given systems |
|
| for(int i = 0; i < List_Nbr(TLAsystems_L); i++){ | | for(int i = 0; i < List_Nbr(TLAsystems_L); i++) { |
| List_Read(TLAsystems_L, i, &TLAsystem); | | List_Read(TLAsystems_L, i, &TLAsystem); |
|
| System_P = (DefineSystem*)List_Pointer(Resolution_P->DefineSystem, | | System_P = (DefineSystem *)List_Pointer(Resolution_P->DefineSystem, |
| TLAsystem.SystemIndex); | | TLAsystem.SystemIndex); |
| DofData_P = DofData_P0 + TLAsystem.SystemIndex; | | DofData_P = DofData_P0 + TLAsystem.SystemIndex; |
| | |
|
| if (!List_Nbr(DofData_P->Solutions)) | | if(!List_Nbr(DofData_P->Solutions)) |
| Message::Error("No initial solution for system %s", System_P->Name); | | Message::Error("No initial solution for system %s", System_P->Name); |
| | |
|
| Solution_P = (struct Solution*) | | Solution_P = (struct Solution *)List_Pointer( |
| List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1); | | DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1); |
| | |
| TimeStep = (int)Current.TimeStep; | | TimeStep = (int)Current.TimeStep; |
|
| if (Solution_P->TimeStep != TimeStep) { // if we compute a new time step | | if(Solution_P->TimeStep != TimeStep) { // if we compute a new time step |
| Solution_S.TimeStep = TimeStep ; | | Solution_S.TimeStep = TimeStep; |
| Solution_S.Time = Current.Time ; | | Solution_S.Time = Current.Time; |
| Solution_S.TimeImag = Current.TimeImag ; | | Solution_S.TimeImag = Current.TimeImag; |
| Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; | | Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P); |
| Solution_S.SolutionExist = 1 ; | | Solution_S.SolutionExist = 1; |
| LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof); | | LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof); |
| | |
| List_Add(DofData_P->Solutions, &Solution_S); | | List_Add(DofData_P->Solutions, &Solution_S); |
|
| DofData_P->CurrentSolution = (struct Solution*) | | DofData_P->CurrentSolution = (struct Solution *)List_Pointer( |
| List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ; | | DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1); |
| Solution_P = DofData_P->CurrentSolution; | | Solution_P = DofData_P->CurrentSolution; |
| } | | } |
| else { | | else { |
| // fix time values if we recompute the same step (with different time) | | // fix time values if we recompute the same step (with different time) |
|
| Solution_P->Time = Current.Time ; | | Solution_P->Time = Current.Time; |
| Solution_P->TimeImag = Current.TimeImag ; | | Solution_P->TimeImag = Current.TimeImag; |
| Free(Solution_P->TimeFunctionValues); | | Free(Solution_P->TimeFunctionValues); |
|
| Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ; | | Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P); |
| } | | } |
| | |
| NbrSolutions = List_Nbr(DofData_P->Solutions); | | NbrSolutions = List_Nbr(DofData_P->Solutions); |
| if(NbrSolutions < Current.PredOrder + 2) | | if(NbrSolutions < Current.PredOrder + 2) |
| Message::Error("Too few past solutions for system %s", System_P->Name); | | Message::Error("Too few past solutions for system %s", System_P->Name); |
| | |
| LinAlg_ZeroVector(&Solution_P->x); | | LinAlg_ZeroVector(&Solution_P->x); |
|
| for (int j=0; j<=Current.PredOrder; j++) { | | for(int j = 0; j <= Current.PredOrder; j++) { |
| PastSolution_P = (struct Solution*)List_Pointer(DofData_P->Solutions, | | PastSolution_P = (struct Solution *)List_Pointer(DofData_P->Solutions, |
| NbrSolutions-2-j); | | NbrSolutions - 2 - j); |
| if (!PastSolution_P->SolutionExist) | | if(!PastSolution_P->SolutionExist) |
| Message::Error("Too few past solutions for system %s", System_P->Name); | | Message::Error("Too few past solutions for system %s", System_P->Name); |
| | |
| x_NminusJ_P = &PastSolution_P->x; | | x_NminusJ_P = &PastSolution_P->x; |
| LinAlg_AddVectorProdVectorDouble(&Solution_P->x, x_NminusJ_P, | | LinAlg_AddVectorProdVectorDouble(&Solution_P->x, x_NminusJ_P, |
| Current.aPredCoeff[j], &Solution_P->x); | | Current.aPredCoeff[j], &Solution_P->x); |
| } | | } |
| | |
|
| xPredicted_P = (gVector*)List_Pointer(xPredicted_L, i); | | xPredicted_P = (gVector *)List_Pointer(xPredicted_L, i); |
| LinAlg_CopyVector(&Solution_P->x, xPredicted_P); | | LinAlg_CopyVector(&Solution_P->x, xPredicted_P); |
| } | | } |
| | |
| // Loop through all specified PostOperations | | // Loop through all specified PostOperations |
|
| if (List_Nbr(TLAPostOp_L) != List_Nbr(Current.PostOpData_L)) | | if(List_Nbr(TLAPostOp_L) != List_Nbr(Current.PostOpData_L)) |
| Message::Error("Current.PostOpData_L list is not up to date"); | | Message::Error("Current.PostOpData_L list is not up to date"); |
|
| for(int i = 0; i < List_Nbr(TLAPostOp_L); i++){ | | for(int i = 0; i < List_Nbr(TLAPostOp_L); i++) { |
| | | PostOpSolutions_P = |
| PostOpSolutions_P = (struct PostOpSolutions*) | | (struct PostOpSolutions *)List_Pointer(Current.PostOpData_L, i); |
| List_Pointer(Current.PostOpData_L, i); | | |
| NbrSolutions = List_Nbr(PostOpSolutions_P->Solutions_L); | | NbrSolutions = List_Nbr(PostOpSolutions_P->Solutions_L); |
|
| if (!NbrSolutions) | | if(!NbrSolutions) |
| Message::Error("No initial result for PostOperation %s", | | Message::Error("No initial result for PostOperation %s", |
| PostOpSolutions_P->PostOperation_P->Name); | | PostOpSolutions_P->PostOperation_P->Name); |
| | |
|
| Solution_P = (struct Solution*)List_Pointer(PostOpSolutions_P->Solutions_L, | | Solution_P = (struct Solution *)List_Pointer(PostOpSolutions_P->Solutions_L, |
| NbrSolutions-1); | | NbrSolutions - 1); |
| | |
| TimeStep = (int)Current.TimeStep; | | TimeStep = (int)Current.TimeStep; |
|
| if (Solution_P->TimeStep != TimeStep) { // if we compute a new time step | | if(Solution_P->TimeStep != TimeStep) { // if we compute a new time step |
| Solution_S.TimeStep = TimeStep ; | | Solution_S.TimeStep = TimeStep; |
| Solution_S.Time = Current.Time ; | | Solution_S.Time = Current.Time; |
| Solution_S.TimeImag = Current.TimeImag ; | | Solution_S.TimeImag = Current.TimeImag; |
| Solution_S.SolutionExist = 1 ; | | Solution_S.SolutionExist = 1; |
| Solution_S.TimeFunctionValues = NULL; | | Solution_S.TimeFunctionValues = NULL; |
| LinAlg_GetVectorSize(&Solution_P->x, &PostOpSolLength); | | LinAlg_GetVectorSize(&Solution_P->x, &PostOpSolLength); |
| LinAlg_CreateVector(&Solution_S.x, &DofData_P0->Solver, PostOpSolLength); | | LinAlg_CreateVector(&Solution_S.x, &DofData_P0->Solver, PostOpSolLength); |
| | |
| List_Add(PostOpSolutions_P->Solutions_L, &Solution_S); | | List_Add(PostOpSolutions_P->Solutions_L, &Solution_S); |
|
| Solution_P = (struct Solution*) | | Solution_P = (struct Solution *)List_Pointer( |
| List_Pointer(PostOpSolutions_P->Solutions_L, | | PostOpSolutions_P->Solutions_L, |
| List_Nbr(PostOpSolutions_P->Solutions_L)-1); | | List_Nbr(PostOpSolutions_P->Solutions_L) - 1); |
| } | | } |
| else { | | else { |
| // fix time values if we recompute the same step (with different time) | | // fix time values if we recompute the same step (with different time) |
|
| Solution_P->Time = Current.Time ; | | Solution_P->Time = Current.Time; |
| Solution_P->TimeImag = Current.TimeImag ; | | Solution_P->TimeImag = Current.TimeImag; |
| } | | } |
| | |
| NbrSolutions = List_Nbr(PostOpSolutions_P->Solutions_L); | | NbrSolutions = List_Nbr(PostOpSolutions_P->Solutions_L); |
| if(NbrSolutions < Current.PredOrder + 2) | | if(NbrSolutions < Current.PredOrder + 2) |
| Message::Error("Too few past results for PostOperation %s", | | Message::Error("Too few past results for PostOperation %s", |
| PostOpSolutions_P->PostOperation_P->Name); | | PostOpSolutions_P->PostOperation_P->Name); |
| | |
|
| PostOpSolPredicted_P = (gVector*)List_Pointer(PostOpSolPredicted_L, i); | | PostOpSolPredicted_P = (gVector *)List_Pointer(PostOpSolPredicted_L, i); |
| LinAlg_ZeroVector(PostOpSolPredicted_P); | | LinAlg_ZeroVector(PostOpSolPredicted_P); |
|
| for (int j=0; j<=Current.PredOrder; j++) { | | for(int j = 0; j <= Current.PredOrder; j++) { |
| PastSolution_P = (struct Solution*)List_Pointer(PostOpSolutions_P->Solutio | | PastSolution_P = (struct Solution *)List_Pointer( |
| ns_L, | | PostOpSolutions_P->Solutions_L, NbrSolutions - 2 - j); |
| NbrSolutions-2-j); | | if(!PastSolution_P->SolutionExist) |
| if (!PastSolution_P->SolutionExist) | | |
| Message::Error("Too few past results for PostOperation %s", | | Message::Error("Too few past results for PostOperation %s", |
|
| PostOpSolutions_P->PostOperation_P->Name); | | PostOpSolutions_P->PostOperation_P->Name); |
| x_NminusJ_P = &PastSolution_P->x; | | x_NminusJ_P = &PastSolution_P->x; |
| LinAlg_AddVectorProdVectorDouble(PostOpSolPredicted_P, x_NminusJ_P, | | LinAlg_AddVectorProdVectorDouble(PostOpSolPredicted_P, x_NminusJ_P, |
|
| Current.aPredCoeff[j], PostOpSolPredicted | | Current.aPredCoeff[j], |
| _P); | | PostOpSolPredicted_P); |
| } | | } |
| } | | } |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l M a x L T E r a t i o */ | | /* C a l M a x L T E r a t i o */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| double CalcMaxLTEratio(Resolution *Resolution_P, | | double CalcMaxLTEratio(Resolution *Resolution_P, DofData *DofData_P0, |
| DofData *DofData_P0, | | List_T *TLAsystems_L, List_T *TLAPostOp_L, int Order, |
| List_T *TLAsystems_L, | | List_T *xPredicted_L, List_T *PostOpSolPredicted_L) |
| List_T *TLAPostOp_L, | | |
| int Order, | | |
| List_T *xPredicted_L, | | |
| List_T *PostOpSolPredicted_L) | | |
| { | | { |
|
| DefineSystem *DefineSystem_P; | | DefineSystem *DefineSystem_P; |
| DofData *DofData_P; | | DofData *DofData_P; |
| PostOpSolutions *PostOpSolutions_P; | | PostOpSolutions *PostOpSolutions_P; |
| TimeLoopAdaptiveSystem TLAsystem; | | TimeLoopAdaptiveSystem TLAsystem; |
| LoopErrorPostOperation TLAPostOp; | | LoopErrorPostOperation TLAPostOp; |
|
| Solution *Solution_P; | | Solution *Solution_P; |
| gVector *xPredictor_P, *xCorrector_P; // predicted and | | gVector *xPredictor_P, *xCorrector_P; // predicted and actual solution vector |
| actual solution vector | | gVector xLTE; // Local Truncation Error vector |
| gVector xLTE; // Local Truncat | | gVector *PostOpSolPred_P, |
| ion Error vector | | *PostOpSolCorr_P; // predicted and actual solution vector |
| gVector *PostOpSolPred_P, *PostOpSolCorr_P; // predicted and | | gVector PostOpSolLTE; // Local Truncation Error vector |
| actual solution vector | | double pec, cec; // predictor and corrector error constants |
| gVector PostOpSolLTE; // Local Truncat | | double ErrorRatio, MaxErrorRatio; |
| ion Error vector | | int NbrSolutions, PostOpSolLength; |
| double pec, cec; // predictor and | | |
| corrector error constants | | |
| double ErrorRatio, MaxErrorRatio; | | |
| int NbrSolutions, PostOpSolLength; | | |
| | |
| MaxErrorRatio = 0.; | | MaxErrorRatio = 0.; |
| | |
| // Determine error constants | | // Determine error constants |
|
| pec = Current.PredErrorConst; // Predictor error constant | | pec = Current.PredErrorConst; // Predictor error constant |
| cec = Current.CorrErrorConst; // Corrector error constant | | cec = Current.CorrErrorConst; // Corrector error constant |
| | |
| // Loop through all given systems | | // Loop through all given systems |
| for(int i = 0; i < List_Nbr(TLAsystems_L); i++) { | | for(int i = 0; i < List_Nbr(TLAsystems_L); i++) { |
| List_Read(TLAsystems_L, i, &TLAsystem); | | List_Read(TLAsystems_L, i, &TLAsystem); |
|
| DefineSystem_P = (DefineSystem*)List_Pointer(Resolution_P->DefineSystem, | | DefineSystem_P = (DefineSystem *)List_Pointer(Resolution_P->DefineSystem, |
| TLAsystem.SystemIndex); | | TLAsystem.SystemIndex); |
| DofData_P = DofData_P0 + TLAsystem.SystemIndex; | | DofData_P = DofData_P0 + TLAsystem.SystemIndex; |
| | |
| NbrSolutions = List_Nbr(DofData_P->Solutions); | | NbrSolutions = List_Nbr(DofData_P->Solutions); |
| if(NbrSolutions < Order + 1) | | if(NbrSolutions < Order + 1) |
|
| Message::Error("Too few past solutions for system %s", DefineSystem_P->Nam | | Message::Error("Too few past solutions for system %s", |
| e); | | DefineSystem_P->Name); |
| | |
|
| xPredictor_P = (gVector*)List_Pointer(xPredicted_L, i); | | xPredictor_P = (gVector *)List_Pointer(xPredicted_L, i); |
| xCorrector_P = &((struct Solution*)List_Pointer(DofData_P->Solutions, | | xCorrector_P = |
| NbrSolutions-1))->x; | | &((struct Solution *)List_Pointer(DofData_P->Solutions, NbrSolutions - 1)) |
| | ->x; |
| | |
| // Vector of all local truncation errors | | // Vector of all local truncation errors |
| // xLTE = cec / (pec - cec) * (xCorrector - xPredictor) | | // xLTE = cec / (pec - cec) * (xCorrector - xPredictor) |
| LinAlg_CreateVector(&xLTE, &DofData_P->Solver, DofData_P->NbrDof); | | LinAlg_CreateVector(&xLTE, &DofData_P->Solver, DofData_P->NbrDof); |
| LinAlg_CopyVector(xCorrector_P, &xLTE); | | LinAlg_CopyVector(xCorrector_P, &xLTE); |
| LinAlg_SubVectorVector(&xLTE, xPredictor_P, &xLTE); | | LinAlg_SubVectorVector(&xLTE, xPredictor_P, &xLTE); |
| LinAlg_ProdVectorDouble(&xLTE, cec / (pec - cec), &xLTE); | | LinAlg_ProdVectorDouble(&xLTE, cec / (pec - cec), &xLTE); |
| | |
|
| Cal_SolutionErrorRatio(&xLTE, xCorrector_P, | | Cal_SolutionErrorRatio(&xLTE, xCorrector_P, TLAsystem.SystemLTEreltol, |
| TLAsystem.SystemLTEreltol, TLAsystem.SystemLTEabstol, | | TLAsystem.SystemLTEabstol, TLAsystem.NormType, |
| TLAsystem.NormType, &ErrorRatio); | | &ErrorRatio); |
| | |
| LinAlg_DestroyVector(&xLTE); | | LinAlg_DestroyVector(&xLTE); |
| | |
|
| if (ErrorRatio != ErrorRatio) { // If ErrorRatio = NaN => There was no v | | if(ErrorRatio != |
| alid solution! | | ErrorRatio) { // If ErrorRatio = NaN => There was no valid solution! |
| MaxErrorRatio = ErrorRatio; | | MaxErrorRatio = ErrorRatio; |
| break; | | break; |
| } | | } |
| | |
|
| if (ErrorRatio > MaxErrorRatio) | | if(ErrorRatio > MaxErrorRatio) MaxErrorRatio = ErrorRatio; |
| MaxErrorRatio = ErrorRatio; | | |
| | |
|
| if(Message::GetVerbosity() > 5) | | if(Message::GetVerbosity() > 5) { |
| { | | |
| Message::Info("LTE %s of error ratio from system %s: %.3g", | | Message::Info("LTE %s of error ratio from system %s: %.3g", |
|
| TLAsystem.NormTypeString, DefineSystem_P->Name, ErrorRatio); | | TLAsystem.NormTypeString, DefineSystem_P->Name, ErrorRatio); |
| } | | } |
| } | | } |
| | |
| // Loop through all given PostOperations | | // Loop through all given PostOperations |
|
| if (List_Nbr(TLAPostOp_L) != List_Nbr(Current.PostOpData_L)) | | if(List_Nbr(TLAPostOp_L) != List_Nbr(Current.PostOpData_L)) |
| Message::Error("Current PostOpData_L list is not up to date"); | | Message::Error("Current PostOpData_L list is not up to date"); |
| for(int i = 0; i < List_Nbr(TLAPostOp_L); i++) { | | for(int i = 0; i < List_Nbr(TLAPostOp_L); i++) { |
| List_Read(TLAPostOp_L, i, &TLAPostOp); | | List_Read(TLAPostOp_L, i, &TLAPostOp); |
|
| PostOpSolutions_P = (struct PostOpSolutions*) | | PostOpSolutions_P = |
| List_Pointer(Current.PostOpData_L, i); | | (struct PostOpSolutions *)List_Pointer(Current.PostOpData_L, i); |
| NbrSolutions = List_Nbr(PostOpSolutions_P->Solutions_L); | | NbrSolutions = List_Nbr(PostOpSolutions_P->Solutions_L); |
| if(NbrSolutions < Order + 1) | | if(NbrSolutions < Order + 1) |
| Message::Error("Too few past solutions for PostOperations %s", | | Message::Error("Too few past solutions for PostOperations %s", |
| PostOpSolutions_P->PostOperation_P->Name); | | PostOpSolutions_P->PostOperation_P->Name); |
|
| Solution_P = (struct Solution*) | | Solution_P = (struct Solution *)List_Pointer(PostOpSolutions_P->Solutions_L, |
| List_Pointer(PostOpSolutions_P->Solutions_L, NbrSolutions-1) | | NbrSolutions - 1); |
| ; | | |
| | |
|
| PostOpSolPred_P = (gVector*)List_Pointer(PostOpSolPredicted_L, i); | | PostOpSolPred_P = (gVector *)List_Pointer(PostOpSolPredicted_L, i); |
| PostOpSolCorr_P = &Solution_P->x; | | PostOpSolCorr_P = &Solution_P->x; |
| | |
| // Vector of all local truncation errors | | // Vector of all local truncation errors |
| // xLTE = cec / (pec - cec) * (xCorrector - xPredictor) | | // xLTE = cec / (pec - cec) * (xCorrector - xPredictor) |
| LinAlg_GetVectorSize(PostOpSolCorr_P, &PostOpSolLength); | | LinAlg_GetVectorSize(PostOpSolCorr_P, &PostOpSolLength); |
| LinAlg_CreateVector(&PostOpSolLTE, &DofData_P0->Solver, PostOpSolLength); | | LinAlg_CreateVector(&PostOpSolLTE, &DofData_P0->Solver, PostOpSolLength); |
| LinAlg_CopyVector(PostOpSolCorr_P, &PostOpSolLTE); | | LinAlg_CopyVector(PostOpSolCorr_P, &PostOpSolLTE); |
| LinAlg_SubVectorVector(&PostOpSolLTE, PostOpSolPred_P, &PostOpSolLTE); | | LinAlg_SubVectorVector(&PostOpSolLTE, PostOpSolPred_P, &PostOpSolLTE); |
| LinAlg_ProdVectorDouble(&PostOpSolLTE, cec / (pec - cec), &PostOpSolLTE); | | LinAlg_ProdVectorDouble(&PostOpSolLTE, cec / (pec - cec), &PostOpSolLTE); |
| | |
|
| Cal_SolutionErrorRatio(&PostOpSolLTE, PostOpSolCorr_P, | | Cal_SolutionErrorRatio( |
| TLAPostOp.PostOperationReltol, TLAPostOp.PostOperatio | | &PostOpSolLTE, PostOpSolCorr_P, TLAPostOp.PostOperationReltol, |
| nAbstol, | | TLAPostOp.PostOperationAbstol, TLAPostOp.NormType, &ErrorRatio); |
| TLAPostOp.NormType, &ErrorRatio); | | |
| | |
| LinAlg_DestroyVector(&PostOpSolLTE); | | LinAlg_DestroyVector(&PostOpSolLTE); |
| | |
|
| if (ErrorRatio != ErrorRatio) { // If ErrorRatio = NaN => There was no v | | if(ErrorRatio != |
| alid solution! | | ErrorRatio) { // If ErrorRatio = NaN => There was no valid solution! |
| MaxErrorRatio = ErrorRatio; | | MaxErrorRatio = ErrorRatio; |
| break; | | break; |
| } | | } |
| | |
|
| if (ErrorRatio > MaxErrorRatio) | | if(ErrorRatio > MaxErrorRatio) MaxErrorRatio = ErrorRatio; |
| MaxErrorRatio = ErrorRatio; | | |
| | |
|
| if(Message::GetVerbosity() > 5) | | if(Message::GetVerbosity() > 5) { |
| { | | |
| Message::Info("LTE %s of error ratio from PostOperation %s: %.3g", | | Message::Info("LTE %s of error ratio from PostOperation %s: %.3g", |
| TLAPostOp.NormTypeString, | | TLAPostOp.NormTypeString, |
| PostOpSolutions_P->PostOperation_P->Name, ErrorRatio); | | PostOpSolutions_P->PostOperation_P->Name, ErrorRatio); |
| } | | } |
| } | | } |
| | |
| return MaxErrorRatio; | | return MaxErrorRatio; |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* G e t I n t e g r a t i o n S c h e m e */ | | /* G e t I n t e g r a t i o n S c h e m e */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void GetIntegrationScheme(Operation *Operation_P, | | void GetIntegrationScheme(Operation *Operation_P, int *TypeTime, int *MaxOrder) |
| int *TypeTime, | | |
| int *MaxOrder) | | |
| { | | { |
|
| if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Euler")) { | | if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Euler")) { |
| *TypeTime = TIME_THETA; | | *TypeTime = TIME_THETA; |
| *MaxOrder = 1; | | *MaxOrder = 1; |
| } | | } |
|
| else if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Trapezoidal")) { | | else if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Trapezoidal")) { |
| *TypeTime = TIME_THETA; | | *TypeTime = TIME_THETA; |
| *MaxOrder = 2; | | *MaxOrder = 2; |
| } | | } |
|
| else if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_2") || | | else if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_2") || |
| !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_2")) { | | !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_2")) { |
| *TypeTime = TIME_GEAR; | | *TypeTime = TIME_GEAR; |
| *MaxOrder = 2; | | *MaxOrder = 2; |
| } | | } |
|
| else if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_3") || | | else if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_3") || |
| !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_3")) { | | !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_3")) { |
| *TypeTime = TIME_GEAR; | | *TypeTime = TIME_GEAR; |
| *MaxOrder = 3; | | *MaxOrder = 3; |
| } | | } |
|
| else if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_4") || | | else if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_4") || |
| !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_4")) { | | !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_4")) { |
| *TypeTime = TIME_GEAR; | | *TypeTime = TIME_GEAR; |
| *MaxOrder = 4; | | *MaxOrder = 4; |
| } | | } |
|
| else if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_5") || | | else if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_5") || |
| !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_5")) { | | !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_5")) { |
| *TypeTime = TIME_GEAR; | | *TypeTime = TIME_GEAR; |
| *MaxOrder = 5; | | *MaxOrder = 5; |
| } | | } |
|
| else if (!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_6") || | | else if(!strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "Gear_6") || |
| !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_6")) { | | !strcmp(Operation_P->Case.TimeLoopAdaptive.Scheme, "BDF_6")) { |
| *TypeTime = TIME_GEAR; | | *TypeTime = TIME_GEAR; |
| *MaxOrder = 6; | | *MaxOrder = 6; |
| } | | } |
| else | | else |
| Message::Error("Unknown integration scheme: %s", | | Message::Error("Unknown integration scheme: %s", |
| Operation_P->Case.TimeLoopAdaptive.Scheme); | | Operation_P->Case.TimeLoopAdaptive.Scheme); |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* O p e r a t i o n _ T i m e L o o p A d a p t i v e */ | | /* O p e r a t i o n _ T i m e L o o p A d a p t i v e */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Operation_TimeLoopAdaptive(Resolution *Resolution_P, | | void Operation_TimeLoopAdaptive(Resolution *Resolution_P, |
| Operation *Operation_P, | | Operation *Operation_P, DofData *DofData_P0, |
| DofData *DofData_P0, | | GeoData *GeoData_P0, int *Flag_Break) |
| GeoData *GeoData_P0, | | |
| int *Flag_Break) | | |
| { | | { |
|
| int TypeTime=0, MaxOrder=0, Order=0, TLATimeStep; | | int TypeTime = 0, MaxOrder = 0, Order = 0, TLATimeStep; |
| int Try, BreakpointNum, NbrSolutions=0, NbrPostOps; | | int Try, BreakpointNum, NbrSolutions = 0, NbrPostOps; |
| | |
|
| double Save_Time, Save_DTime, Save_Theta, maxLTEratio=0, nextBreakpoint; | | double Save_Time, Save_DTime, Save_Theta, maxLTEratio = 0, nextBreakpoint; |
| double Save_TimeStep, FirstTimePoint, DTimeBeforeBreakpoint=1.; | | double Save_TimeStep, FirstTimePoint, DTimeBeforeBreakpoint = 1.; |
| bool TimeStepAccepted=true, DTimeMinAtLastStep, BreakpointListCreated; | | bool TimeStepAccepted = true, DTimeMinAtLastStep, BreakpointListCreated; |
| bool BreakpointAtThisStep, BreakpointAtNextStep; | | bool BreakpointAtThisStep, BreakpointAtNextStep; |
| double Time0, TimeMax, DTimeInit, DTimeMin, DTimeMax; | | double Time0, TimeMax, DTimeInit, DTimeMin, DTimeMax; |
| double LTEtarget, DTimeMaxScal, DTimeScal_NotConverged, DTimeScal_PETScError; | | double LTEtarget, DTimeMaxScal, DTimeScal_NotConverged, DTimeScal_PETScError; |
|
| double DTimeScal=1.0; | | double DTimeScal = 1.0; |
| List_T *Breakpoints_L, *TLAsystems_L, *LEPostOp_L; | | List_T *Breakpoints_L, *TLAsystems_L, *LEPostOp_L; |
| List_T *LEPostOpNames_L; | | List_T *LEPostOpNames_L; |
| List_T *xPredicted_L, *PostOpSolPredicted_L; | | List_T *xPredicted_L, *PostOpSolPredicted_L; |
| TimeLoopAdaptiveSystem TLAsystem; | | TimeLoopAdaptiveSystem TLAsystem; |
|
| DofData *DofData_P=NULL; | | DofData *DofData_P = NULL; |
| gVector xPredicted_S; | | gVector xPredicted_S; |
| | |
| // Some default values for constants influencing the time stepping | | // Some default values for constants influencing the time stepping |
|
| LTEtarget = 0.8; // target LTE ratio for next step (should be | | LTEtarget = 0.8; // target LTE ratio for next step (should be below 1) |
| below 1) | | DTimeMaxScal = 2.0; // maximum factor for increasing the time step DTime |
| DTimeMaxScal = 2.0; // maximum factor for increasing the time ste | | DTimeScal_NotConverged = |
| p DTime | | 0.25; // step size scaling in case of a not converged iterative loop |
| DTimeScal_NotConverged = 0.25; // step size scaling in case of a not converg | | DTimeScal_PETScError = 0.25; // step size scaling in case of a PETSc error |
| ed iterative loop | | |
| DTimeScal_PETScError = 0.25; // step size scaling in case of a PETSc error | | |
| | |
| // Override default values if they are provided by the user | | // Override default values if they are provided by the user |
| LTEtarget = (Operation_P->Case.TimeLoopAdaptive.LTEtarget < 0) ? | | LTEtarget = (Operation_P->Case.TimeLoopAdaptive.LTEtarget < 0) ? |
|
| LTEtarget : Operation_P->Case.TimeLoopAdaptive.LTEtarget; | | LTEtarget : |
| DTimeMaxScal = (Operation_P->Case.TimeLoopAdaptive.DTimeMaxScal < 0 ) ? | | Operation_P->Case.TimeLoopAdaptive.LTEtarget; |
| DTimeMaxScal : Operation_P->Case.TimeLoopAdaptive.DTimeMaxScal; | | DTimeMaxScal = (Operation_P->Case.TimeLoopAdaptive.DTimeMaxScal < 0) ? |
| | DTimeMaxScal : |
| | Operation_P->Case.TimeLoopAdaptive.DTimeMaxScal; |
| DTimeScal_NotConverged = | | DTimeScal_NotConverged = |
|
| (Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged < 0) ? | | (Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged < 0) ? |
| DTimeScal_NotConverged : | | DTimeScal_NotConverged : |
| Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged; | | Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged; |
| DTimeScal_PETScError = | | DTimeScal_PETScError = |
|
| (Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged < 0) ? | | (Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged < 0) ? |
| DTimeScal_PETScError : | | DTimeScal_PETScError : |
| Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged; | | Operation_P->Case.TimeLoopAdaptive.DTimeScal_NotConverged; |
| | |
|
| Time0 = Operation_P->Case.TimeLoopAdaptive.Time0; | | Time0 = Operation_P->Case.TimeLoopAdaptive.Time0; |
| TimeMax = Operation_P->Case.TimeLoopAdaptive.TimeMax; | | TimeMax = Operation_P->Case.TimeLoopAdaptive.TimeMax; |
| DTimeInit = Operation_P->Case.TimeLoopAdaptive.DTimeInit; | | DTimeInit = Operation_P->Case.TimeLoopAdaptive.DTimeInit; |
| DTimeMin = Operation_P->Case.TimeLoopAdaptive.DTimeMin; | | DTimeMin = Operation_P->Case.TimeLoopAdaptive.DTimeMin; |
| DTimeMax = Operation_P->Case.TimeLoopAdaptive.DTimeMax; | | DTimeMax = Operation_P->Case.TimeLoopAdaptive.DTimeMax; |
| Breakpoints_L = Operation_P->Case.TimeLoopAdaptive.Breakpoints_L; | | Breakpoints_L = Operation_P->Case.TimeLoopAdaptive.Breakpoints_L; |
| TLAsystems_L = Operation_P->Case.TimeLoopAdaptive.TimeLoopAdaptiveSystems_L; | | TLAsystems_L = Operation_P->Case.TimeLoopAdaptive.TimeLoopAdaptiveSystems_L; |
| LEPostOp_L = Operation_P->Case.TimeLoopAdaptive.TimeLoopAdaptivePOs_L; | | LEPostOp_L = Operation_P->Case.TimeLoopAdaptive.TimeLoopAdaptivePOs_L; |
| GetIntegrationScheme(Operation_P, &TypeTime, &MaxOrder); | | GetIntegrationScheme(Operation_P, &TypeTime, &MaxOrder); |
| | |
|
| xPredicted_L = List_Create(4,4,sizeof(gVector)); | | xPredicted_L = List_Create(4, 4, sizeof(gVector)); |
| PostOpSolPredicted_L = List_Create(4,4,sizeof(gVector)); | | PostOpSolPredicted_L = List_Create(4, 4, sizeof(gVector)); |
| | |
| // Just some checks | | // Just some checks |
| // ---------------- | | // ---------------- |
| | |
|
| if (TLAsystems_L == NULL) | | if(TLAsystems_L == NULL) |
| TLAsystems_L = List_Create(1,1,sizeof(TimeLoopAdaptiveSystem)); | | TLAsystems_L = List_Create(1, 1, sizeof(TimeLoopAdaptiveSystem)); |
| if (LEPostOp_L == NULL) | | if(LEPostOp_L == NULL) |
| LEPostOp_L = List_Create(1,1,sizeof(LoopErrorPostOperation)); | | LEPostOp_L = List_Create(1, 1, sizeof(LoopErrorPostOperation)); |
| | |
| // Check the timing values | | // Check the timing values |
|
| if (Time0 > TimeMax) | | if(Time0 > TimeMax) Message::Error("Time0 > TimeMax"); |
| Message::Error("Time0 > TimeMax"); | | if(DTimeInit < DTimeMin) Message::Error("DTimeInit < DTimeMin"); |
| if (DTimeInit < DTimeMin) | | if(DTimeInit > DTimeMax) Message::Error("DTimeInit > DTimeMax"); |
| Message::Error("DTimeInit < DTimeMin"); | | if(DTimeInit > TimeMax - Time0) |
| if (DTimeInit > DTimeMax) | | Message::Error("DTimeInit > (TimeMax - Time0"); |
| Message::Error("DTimeInit > DTimeMax"); | | |
| if (DTimeInit > TimeMax - Time0) | | |
| Message::Error("DTimeInit > (TimeMax - Time0"); | | |
| | |
| // Initialization before starting the time loop | | // Initialization before starting the time loop |
| // -------------------------------------------- | | // -------------------------------------------- |
| | |
| // Check if initial solutions for all specified systems are available | | // Check if initial solutions for all specified systems are available |
| // and create vectors for the predicted solutions | | // and create vectors for the predicted solutions |
|
| for(int i = 0; i < List_Nbr(TLAsystems_L); i++){ | | for(int i = 0; i < List_Nbr(TLAsystems_L); i++) { |
| List_Read(TLAsystems_L, i, &TLAsystem); | | List_Read(TLAsystems_L, i, &TLAsystem); |
|
| DefineSystem *System_P = (DefineSystem*)List_Pointer(Resolution_P->DefineSys | | DefineSystem *System_P = (DefineSystem *)List_Pointer( |
| tem, | | Resolution_P->DefineSystem, TLAsystem.SystemIndex); |
| TLAsystem.SystemIndex); | | |
| DofData_P = DofData_P0 + TLAsystem.SystemIndex; | | DofData_P = DofData_P0 + TLAsystem.SystemIndex; |
| NbrSolutions = List_Nbr(DofData_P->Solutions); | | NbrSolutions = List_Nbr(DofData_P->Solutions); |
| | |
| if(!NbrSolutions) | | if(!NbrSolutions) |
| Message::Error("No initial solution for system %s", System_P->Name); | | Message::Error("No initial solution for system %s", System_P->Name); |
| | |
| LinAlg_CreateVector(&xPredicted_S, &DofData_P->Solver, DofData_P->NbrDof); | | LinAlg_CreateVector(&xPredicted_S, &DofData_P->Solver, DofData_P->NbrDof); |
| List_Add(xPredicted_L, &xPredicted_S); | | List_Add(xPredicted_L, &xPredicted_S); |
| } | | } |
| | |
| // Initializing stuff for PostOperations | | // Initializing stuff for PostOperations |
| NbrPostOps = List_Nbr(LEPostOp_L); | | NbrPostOps = List_Nbr(LEPostOp_L); |
|
| LEPostOpNames_L = List_Create(NbrPostOps,1,sizeof(char *)); | | LEPostOpNames_L = List_Create(NbrPostOps, 1, sizeof(char *)); |
| InitLEPostOperation(Resolution_P, DofData_P0, GeoData_P0, LEPostOp_L, | | InitLEPostOperation(Resolution_P, DofData_P0, GeoData_P0, LEPostOp_L, |
| LEPostOpNames_L, PostOpSolPredicted_L); | | LEPostOpNames_L, PostOpSolPredicted_L); |
| | |
| // Some other necessary initializations | | // Some other necessary initializations |
| if(Flag_RESTART && NbrSolutions > 1) | | if(Flag_RESTART && NbrSolutions > 1) |
|
| Current.DTime = ((struct Solution*)List_Pointer(DofData_P->Solutions, | | Current.DTime = |
| NbrSolutions-1))->Time - | | ((struct Solution *)List_Pointer(DofData_P->Solutions, NbrSolutions - 1)) |
| ((struct Solution*)List_Pointer(DofData_P->Solutions, | | ->Time - |
| NbrSolutions-2))->Time; | | ((struct Solution *)List_Pointer(DofData_P->Solutions, NbrSolutions - 2)) |
| | ->Time; |
| else | | else |
| Current.DTime = DTimeInit; | | Current.DTime = DTimeInit; |
| | |
| if(Flag_RESTART) { | | if(Flag_RESTART) { |
|
| if (Current.Time < TimeMax) | | if(Current.Time < TimeMax) Flag_RESTART = 0; |
| Flag_RESTART = 0 ; | | |
| } | | } |
| else | | else |
|
| Current.Time = Time0 ; | | Current.Time = Time0; |
| | |
| Current.TimeStep += 1.0; | | Current.TimeStep += 1.0; |
| TLATimeStep = 1; | | TLATimeStep = 1; |
| // Starting with 1st order (Backward Euler corrector) | | // Starting with 1st order (Backward Euler corrector) |
| Order = 1; | | Order = 1; |
|
| if (TypeTime == TIME_THETA) | | if(TypeTime == TIME_THETA) Current.Theta = 1.0; |
| Current.Theta = 1.0; | | |
| | |
| BreakpointListCreated = !Breakpoints_L; | | BreakpointListCreated = !Breakpoints_L; |
|
| if (BreakpointListCreated) | | if(BreakpointListCreated) Breakpoints_L = List_Create(1, 1, sizeof(double)); |
| Breakpoints_L = List_Create(1,1,sizeof(double)); | | |
| List_Add(Breakpoints_L, &TimeMax); | | List_Add(Breakpoints_L, &TimeMax); |
| List_Sort(Breakpoints_L, fcmp_double); | | List_Sort(Breakpoints_L, fcmp_double); |
| BreakpointNum = 0; | | BreakpointNum = 0; |
| BreakpointAtNextStep = false; | | BreakpointAtNextStep = false; |
| List_Read(Breakpoints_L, BreakpointNum, &nextBreakpoint); | | List_Read(Breakpoints_L, BreakpointNum, &nextBreakpoint); |
|
| FirstTimePoint = Current.Time+Current.DTime; | | FirstTimePoint = Current.Time + Current.DTime; |
| Current.Breakpoint = List_ISearchSeq(Breakpoints_L, &FirstTimePoint, fcmp_doub | | Current.Breakpoint = |
| le); | | List_ISearchSeq(Breakpoints_L, &FirstTimePoint, fcmp_double); |
| if (Current.Breakpoint >= 0) { | | if(Current.Breakpoint >= 0) { |
| BreakpointAtNextStep = true; | | BreakpointAtNextStep = true; |
| DTimeBeforeBreakpoint = Current.DTime; | | DTimeBeforeBreakpoint = Current.DTime; |
| } | | } |
|
| for (int i = 0; i < List_Nbr(Breakpoints_L); i++){ | | for(int i = 0; i < List_Nbr(Breakpoints_L); i++) { |
| List_Read(Breakpoints_L, i, &nextBreakpoint); | | List_Read(Breakpoints_L, i, &nextBreakpoint); |
|
| if (nextBreakpoint > (FirstTimePoint + DTimeMin)) { | | if(nextBreakpoint > (FirstTimePoint + DTimeMin)) { |
| BreakpointNum = i; | | BreakpointNum = i; |
| break; | | break; |
| } | | } |
| } | | } |
| Try = 0; | | Try = 0; |
| | |
| // Start the time loop | | // Start the time loop |
| // ------------------- | | // ------------------- |
| | |
|
| while (Current.Time < TimeMax) { | | while(Current.Time < TimeMax) { |
| | | |
| if(Message::GetOnelabAction() == "stop") break; | | if(Message::GetOnelabAction() == "stop") break; |
| | |
| Message::SetOperatingInTimeLoopAdaptive(true); | | Message::SetOperatingInTimeLoopAdaptive(true); |
| | |
|
| Current.TypeTime = TypeTime; | | Current.TypeTime = TypeTime; |
| Current.Time += Current.DTime; | | Current.Time += Current.DTime; |
| Save_DTime = Current.DTime; | | Save_DTime = Current.DTime; |
| Save_Time = Current.Time; | | Save_Time = Current.Time; |
| Save_TimeStep = Current.TimeStep; | | Save_TimeStep = Current.TimeStep; |
| Save_Theta = Current.Theta; | | Save_Theta = Current.Theta; |
| Try++; | | Try++; |
| BreakpointAtThisStep = BreakpointAtNextStep; | | BreakpointAtThisStep = BreakpointAtNextStep; |
| | |
| Message::SetLastPETScError(0); | | Message::SetLastPETScError(0); |
| | |
|
| Message::Info("Time step %d Try %d Time = %.8g s Stepsize = %.8g s Integ | | Message::Info("Time step %d Try %d Time = %.8g s Stepsize = %.8g s " |
| r. Order = %d", | | "Integr. Order = %d", |
| (int)Current.TimeStep, Try, Current.Time, Current.DTime, Order | | (int)Current.TimeStep, Try, Current.Time, Current.DTime, |
| ); | | Order); |
| if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterStep() < | | if(Message::GetProgressMeterStep() > 0 && |
| 100){ | | Message::GetProgressMeterStep() < 100) { |
| Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + | | Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + |
|
| "/TimeLoopAdaptive/Time", | | "/TimeLoopAdaptive/Time", |
| std::vector<double>(1, Current.Time)); | | std::vector<double>(1, Current.Time)); |
| Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + | | Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + |
|
| "/TimeLoopAdaptive/DTime", | | "/TimeLoopAdaptive/DTime", |
| std::vector<double>(1, Current.DTime)); | | std::vector<double>(1, Current.DTime)); |
| } | | } |
| | |
| // Calculate integration coefficients | | // Calculate integration coefficients |
|
| CalcIntegrationCoefficients(Resolution_P, DofData_P0,TLAsystems_L, | | CalcIntegrationCoefficients(Resolution_P, DofData_P0, TLAsystems_L, |
| LEPostOp_L, Order); | | LEPostOp_L, Order); |
| | |
| // Execute predictor | | // Execute predictor |
| Predictor(Resolution_P, DofData_P0, TLAsystems_L, LEPostOp_L, Order, | | Predictor(Resolution_P, DofData_P0, TLAsystems_L, LEPostOp_L, Order, |
| xPredicted_L, PostOpSolPredicted_L); | | xPredicted_L, PostOpSolPredicted_L); |
| | |
|
| if (NbrPostOps && TimeStepAccepted) | | if(NbrPostOps && TimeStepAccepted) Free_UnusedPOresults(); |
| Free_UnusedPOresults(); | | |
| | |
| // Execute corrector | | // Execute corrector |
| // ----------------- | | // ----------------- |
| Flag_IterativeLoopConverged = 1; | | Flag_IterativeLoopConverged = 1; |
| | |
| Treatment_Operation(Resolution_P, | | Treatment_Operation(Resolution_P, |
| Operation_P->Case.TimeLoopAdaptive.Operation, | | Operation_P->Case.TimeLoopAdaptive.Operation, |
|
| DofData_P0, GeoData_P0, NULL, NULL) ; | | DofData_P0, GeoData_P0, NULL, NULL); |
| Current.Time = Save_Time; | | Current.Time = Save_Time; |
| Current.TypeTime = TypeTime; | | Current.TypeTime = TypeTime; |
|
| Current.DTime = Save_DTime; | | Current.DTime = Save_DTime; |
| Current.TimeStep = Save_TimeStep; | | Current.TimeStep = Save_TimeStep; |
|
| Current.Theta = Save_Theta; | | Current.Theta = Save_Theta; |
| | |
| if(*Flag_Break) { | | if(*Flag_Break) { |
| *Flag_Break = 0; | | *Flag_Break = 0; |
| Message::Info("Flag Break detected. Aborting TimeLoopAdaptive"); | | Message::Info("Flag Break detected. Aborting TimeLoopAdaptive"); |
| break; | | break; |
| } | | } |
| | |
| // Assessing the current time step and eventually | | // Assessing the current time step and eventually |
| // execute the 2nd set of operations | | // execute the 2nd set of operations |
| // ---------------------------------------------- | | // ---------------------------------------------- |
|
| if (Flag_IterativeLoopConverged != 1){ | | if(Flag_IterativeLoopConverged != 1) { |
| TimeStepAccepted = false; | | TimeStepAccepted = false; |
| DTimeScal = DTimeScal_NotConverged; | | DTimeScal = DTimeScal_NotConverged; |
|
| Message::Info("Time step %d Try %d Time = %.8g s rejected (IterativeLoo | | Message::Info( |
| p not " | | "Time step %d Try %d Time = %.8g s rejected (IterativeLoop not " |
| "converged)", (int)Current.TimeStep, Try, Current.Time); | | "converged)", |
| | (int)Current.TimeStep, Try, Current.Time); |
| } | | } |
|
| else if (Message::GetLastPETScError()) { | | else if(Message::GetLastPETScError()) { |
| TimeStepAccepted = false; | | TimeStepAccepted = false; |
| Flag_IterativeLoopConverged = 0; | | Flag_IterativeLoopConverged = 0; |
| DTimeScal = DTimeScal_PETScError; | | DTimeScal = DTimeScal_PETScError; |
| Message::Warning("Time step %d Try %d Time = %.8g s rejected:", | | Message::Warning("Time step %d Try %d Time = %.8g s rejected:", |
| (int)Current.TimeStep, Try, Current.Time); | | (int)Current.TimeStep, Try, Current.Time); |
| Message::Warning("No valid solution found (PETSc-Error: %d)!", | | Message::Warning("No valid solution found (PETSc-Error: %d)!", |
| Message::GetLastPETScError()); | | Message::GetLastPETScError()); |
| Message::SetLastPETScError(0); | | Message::SetLastPETScError(0); |
| } | | } |
|
| else{ | | else { |
| if (NbrPostOps) // Execute the PostOperations if necessary | | if(NbrPostOps) // Execute the PostOperations if necessary |
| Operation_PostOperation(Resolution_P, DofData_P0, GeoData_P0, LEPostOpNa | | Operation_PostOperation(Resolution_P, DofData_P0, GeoData_P0, |
| mes_L); | | LEPostOpNames_L); |
| maxLTEratio = CalcMaxLTEratio(Resolution_P, DofData_P0, TLAsystems_L, LEPo | | maxLTEratio = |
| stOp_L, | | CalcMaxLTEratio(Resolution_P, DofData_P0, TLAsystems_L, LEPostOp_L, |
| Order, xPredicted_L, PostOpSolPredicted_L); | | Order, xPredicted_L, PostOpSolPredicted_L); |
| if (maxLTEratio != maxLTEratio) { // If maxLTEratio = NaN => There was no | | if(maxLTEratio != |
| valid solution! | | maxLTEratio) { // If maxLTEratio = NaN => There was no valid solution! |
| TimeStepAccepted = false; | | TimeStepAccepted = false; |
| Flag_IterativeLoopConverged = 0; | | Flag_IterativeLoopConverged = 0; |
| DTimeScal = DTimeScal_PETScError; | | DTimeScal = DTimeScal_PETScError; |
|
| Message::Info("Time step %d Try %d Time = %.8g s rejected: No valid s | | Message::Info( |
| olution " | | "Time step %d Try %d Time = %.8g s rejected: No valid solution " |
| "found (NaN or Inf)!", (int)Current.TimeStep, Try, Current | | "found (NaN or Inf)!", |
| .Time); | | (int)Current.TimeStep, Try, Current.Time); |
| } | | } |
| else { | | else { |
| if(Message::GetVerbosity() > 4) | | if(Message::GetVerbosity() > 4) |
| Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + | | Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + |
|
| "/TimeLoopAdaptive/LTEmaxErrorRatio", | | "/TimeLoopAdaptive/LTEmaxErrorRatio", |
| std::vector<double>(1, maxLTEratio)); | | std::vector<double>(1, maxLTEratio)); |
|
| if (maxLTEratio <= 1.0){ | | if(maxLTEratio <= 1.0) { |
| TimeStepAccepted = true; | | TimeStepAccepted = true; |
|
| Message::Info("Time step %d Try %d Time = %.8g s accepted (max. LTE | | Message::Info("Time step %d Try %d Time = %.8g s accepted (max. " |
| ratio = %.3g)", | | "LTE ratio = %.3g)", |
| (int)Current.TimeStep, Try, Current.Time, maxLTEratio); | | (int)Current.TimeStep, Try, Current.Time, maxLTEratio); |
| } | | } |
|
| else{ | | else { |
| TimeStepAccepted = false; | | TimeStepAccepted = false; |
|
| Message::Info("Time step %d Try %d Time = %.8g s rejected (max. LTE | | Message::Info("Time step %d Try %d Time = %.8g s rejected (max. " |
| ratio = %.3g)", | | "LTE ratio = %.3g)", |
| (int)Current.TimeStep, Try, Current.Time, maxLTEratio); | | (int)Current.TimeStep, Try, Current.Time, maxLTEratio); |
| } | | } |
| } | | } |
| } | | } |
| | |
|
| if (TimeStepAccepted == true) { | | if(TimeStepAccepted == true) { |
| Treatment_Operation(Resolution_P, | | Treatment_Operation(Resolution_P, |
| Operation_P->Case.TimeLoopAdaptive.OperationEnd, | | Operation_P->Case.TimeLoopAdaptive.OperationEnd, |
|
| DofData_P0, GeoData_P0, NULL, NULL) ; | | DofData_P0, GeoData_P0, NULL, NULL); |
| Current.Time = Save_Time; | | Current.Time = Save_Time; |
| Current.TypeTime = TypeTime; | | Current.TypeTime = TypeTime; |
|
| Current.DTime = Save_DTime; | | Current.DTime = Save_DTime; |
| Current.TimeStep = Save_TimeStep; | | Current.TimeStep = Save_TimeStep; |
|
| Current.Theta = Save_Theta; | | Current.Theta = Save_Theta; |
| Current.TimeStep += 1.; | | Current.TimeStep += 1.; |
|
| TLATimeStep += 1; | | TLATimeStep += 1; |
| Try = 0; | | Try = 0; |
| } | | } |
|
| else{ | | else { |
| if (BreakpointAtThisStep) { | | if(BreakpointAtThisStep) { |
| BreakpointNum = List_ISearchSeq(Breakpoints_L, &Current.Time, fcmp_doubl | | BreakpointNum = |
| e); | | List_ISearchSeq(Breakpoints_L, &Current.Time, fcmp_double); |
| List_Read(Breakpoints_L, BreakpointNum, &nextBreakpoint); | | List_Read(Breakpoints_L, BreakpointNum, &nextBreakpoint); |
| } | | } |
| Current.Time -= Current.DTime; | | Current.Time -= Current.DTime; |
|
| BreakpointAtThisStep = (bool) List_Search(Breakpoints_L, &Current.Time, fc | | BreakpointAtThisStep = |
| mp_double); | | (bool)List_Search(Breakpoints_L, &Current.Time, fcmp_double); |
| } | | } |
| | |
| if(*Flag_Break) { | | if(*Flag_Break) { |
| *Flag_Break = 0; | | *Flag_Break = 0; |
| Message::Info("Flag Break detected. Aborting TimeLoopAdaptive"); | | Message::Info("Flag Break detected. Aborting TimeLoopAdaptive"); |
| break; | | break; |
| } | | } |
| | |
| // Calculate new time step | | // Calculate new time step |
| // ----------------------- | | // ----------------------- |
| DTimeMinAtLastStep = Current.DTime <= DTimeMin; | | DTimeMinAtLastStep = Current.DTime <= DTimeMin; |
|
| if (TimeStepAccepted == false && | | if(TimeStepAccepted == false && DTimeMinAtLastStep && Order < 2) |
| DTimeMinAtLastStep && | | |
| Order < 2) | | |
| Message::Error("Time step too small! Simulation aborted!"); | | Message::Error("Time step too small! Simulation aborted!"); |
| | |
|
| if (Flag_IterativeLoopConverged == 1){ | | if(Flag_IterativeLoopConverged == 1) { |
| // Milne's estimate | | // Milne's estimate |
|
| if (maxLTEratio <= 0) | | if(maxLTEratio <= 0) |
| DTimeScal = DTimeMaxScal; | | DTimeScal = DTimeMaxScal; |
| else { | | else { |
|
| if (Current.TimeStep < 1.5 || (NbrPostOps > 0 && TLATimeStep <= 2) ) | | if(Current.TimeStep < 1.5 || (NbrPostOps > 0 && TLATimeStep <= 2)) |
| // linear adjustment because predictor is of order 0 | | // linear adjustment because predictor is of order 0 |
|
| DTimeScal = LTEtarget/maxLTEratio; | | DTimeScal = LTEtarget / maxLTEratio; |
| else | | else |
|
| DTimeScal = pow(LTEtarget/maxLTEratio, 1./(Order+1.)); | | DTimeScal = pow(LTEtarget / maxLTEratio, 1. / (Order + 1.)); |
| } | | } |
|
| if (DTimeScal >= DTimeMaxScal) { | | if(DTimeScal >= DTimeMaxScal) { |
| if (BreakpointAtThisStep) { | | if(BreakpointAtThisStep) { |
| double dt1, dt2, dtmax; | | double dt1, dt2, dtmax; |
| dt1 = Current.DTime * DTimeMaxScal; | | dt1 = Current.DTime * DTimeMaxScal; |
| dt2 = DTimeBeforeBreakpoint; | | dt2 = DTimeBeforeBreakpoint; |
| dtmax = (dt1 > dt2) ? dt1 : dt2; | | dtmax = (dt1 > dt2) ? dt1 : dt2; |
| DTimeScal = dtmax / Current.DTime; | | DTimeScal = dtmax / Current.DTime; |
| } | | } |
| else | | else |
| DTimeScal = DTimeMaxScal; | | DTimeScal = DTimeMaxScal; |
| } | | } |
| } | | } |
| Current.DTime *= DTimeScal; | | Current.DTime *= DTimeScal; |
| | |
| // Limit the max step size | | // Limit the max step size |
|
| if (Current.DTime > DTimeMax) | | if(Current.DTime > DTimeMax) Current.DTime = DTimeMax; |
| Current.DTime = DTimeMax; | | |
| | |
| // Check that we do not jump over a breakpoint | | // Check that we do not jump over a breakpoint |
|
| if ((Current.DTime + Current.Time >= nextBreakpoint - DTimeMin) && | | if((Current.DTime + Current.Time >= nextBreakpoint - DTimeMin) && |
| (BreakpointNum >= 0)){ | | (BreakpointNum >= 0)) { |
| DTimeBeforeBreakpoint = Current.DTime; | | DTimeBeforeBreakpoint = Current.DTime; |
| Current.DTime = nextBreakpoint - Current.Time; | | Current.DTime = nextBreakpoint - Current.Time; |
| BreakpointAtNextStep = true; | | BreakpointAtNextStep = true; |
| Current.Breakpoint = BreakpointNum; | | Current.Breakpoint = BreakpointNum; |
|
| if (BreakpointNum < List_Nbr(Breakpoints_L)-1){ | | if(BreakpointNum < List_Nbr(Breakpoints_L) - 1) { |
| // There are further breakpoints | | // There are further breakpoints |
| BreakpointNum++; | | BreakpointNum++; |
| List_Read(Breakpoints_L, BreakpointNum, &nextBreakpoint); | | List_Read(Breakpoints_L, BreakpointNum, &nextBreakpoint); |
| } | | } |
| else | | else |
|
| //No further breakpoint | | // No further breakpoint |
| BreakpointNum = -1; | | BreakpointNum = -1; |
| } | | } |
| else { | | else { |
| BreakpointAtNextStep = false; | | BreakpointAtNextStep = false; |
| Current.Breakpoint = -1.; | | Current.Breakpoint = -1.; |
| } | | } |
| | |
| // Limit the min step size | | // Limit the min step size |
|
| if (Current.DTime < DTimeMin) | | if(Current.DTime < DTimeMin) Current.DTime = DTimeMin; |
| Current.DTime = DTimeMin; | | |
| | |
| // Adjust order | | // Adjust order |
| // ------------ | | // ------------ |
|
| if ( Flag_IterativeLoopConverged != 1 || | | if(Flag_IterativeLoopConverged != 1 || |
| // BreakpointAtThisStep || | | // BreakpointAtThisStep || |
| DTimeMinAtLastStep ) | | DTimeMinAtLastStep) |
| Order = 1; | | Order = 1; |
|
| else if ( TLATimeStep > 2 && | | else if(TLATimeStep > 2 && Current.DTime > DTimeMin && TimeStepAccepted && |
| Current.DTime > DTimeMin && | | !BreakpointAtThisStep && Order < MaxOrder) |
| TimeStepAccepted && | | |
| !BreakpointAtThisStep && | | |
| Order < MaxOrder ) | | |
| Order++; | | Order++; |
| | |
|
| if (TypeTime == TIME_THETA) | | if(TypeTime == TIME_THETA) switch(Order) { |
| switch(Order){ | | case 1: |
| case 1: | | Current.Theta = 1.0; // Corrector: Backward Euler |
| Current.Theta = 1.0; // Corrector: Backward Euler | | break; |
| break; | | case 2: |
| case 2: | | Current.Theta = 0.5; // Corrector: Trapezoidal Method |
| Current.Theta = 0.5; // Corrector: Trapezoidal Method | | break; |
| break; | | default: |
| default : | | Message::Error("Order %d not allowed for Theta scheme.", Order); |
| Message::Error("Order %d not allowed for Theta scheme.", Order); | | break; |
| break; | | |
| } | | } |
| } // while loop | | } // while loop |
| | |
| Message::SetOperatingInTimeLoopAdaptive(false); | | Message::SetOperatingInTimeLoopAdaptive(false); |
| | |
|
| Current.TimeStep -= 1.; // Correct the time step counter | | Current.TimeStep -= 1.; // Correct the time step counter |
| | |
| // Finally clean up, destroy vectors and delete lists | | // Finally clean up, destroy vectors and delete lists |
| // -------------------------------------------------- | | // -------------------------------------------------- |
| | |
| for(int i = 0; i < List_Nbr(TLAsystems_L); i++) | | for(int i = 0; i < List_Nbr(TLAsystems_L); i++) |
|
| LinAlg_DestroyVector((gVector*)List_Pointer(xPredicted_L, i)); | | LinAlg_DestroyVector((gVector *)List_Pointer(xPredicted_L, i)); |
| List_Delete(TLAsystems_L); | | List_Delete(TLAsystems_L); |
| List_Delete(xPredicted_L); | | List_Delete(xPredicted_L); |
| | |
| ClearLEPostOperation(Resolution_P, DofData_P0, GeoData_P0, LEPostOp_L, | | ClearLEPostOperation(Resolution_P, DofData_P0, GeoData_P0, LEPostOp_L, |
|
| LEPostOpNames_L, PostOpSolPredicted_L, true); | | LEPostOpNames_L, PostOpSolPredicted_L, true); |
| | |
|
| if (BreakpointListCreated) | | if(BreakpointListCreated) List_Delete(Breakpoints_L); |
| List_Delete(Breakpoints_L); | | |
| } | | } |
| | | |
| End of changes. 154 change blocks. |
|---|
| 436 lines changed or deleted | | 386 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 