Commit 37d62c4c by Kevin Jacques

Merge branch 'RunTimeNLIterVar' into 'master'

RunTimeNLIterVar

See merge request !33
parents 52966f10 fcecb101
Pipeline #1127 passed with stage
in 9 minutes 35 seconds
......@@ -1102,10 +1102,12 @@ struct StringXFunction2Nbr F_Function[] = { /* #Par #Arg */
//{"Update_Jk" , (CAST)F_Update_Jk , -1, 6 }, //kj+++ // NOT USED FOR NOW (26/06/2016)
//{"Update_Jk_sd" , (CAST)F_Update_Jk_sd , -1, 6 }, //kj+++ // NOT USED FOR NOW (26/06/2016)
{"Update_Cell_K" , (CAST)F_Update_Cell_K , -1, 4 }, //kj+++
{"b_Vinch_K" , (CAST)F_b_Vinch_K , -1, 7 }, // 1+3*2=7 //kj+++
{"h_Vinch_K" , (CAST)F_h_Vinch_K , -1, 9 }, // parameter is dimension {2},{3}, 3+3*2=9 //kj+++
{"dbdh_Vinch_K" , (CAST)F_dbdh_Vinch_K , -1, 7 }, // parameter is dimension {2},{3}, 1+3*2=7 //kj+++
{"dhdb_Vinch_K" , (CAST)F_dhdb_Vinch_K , -1, 7 }, // parameter is dimension {2},{3}, 1+3*2=7 //kj+++
{"b_Vinch_K" , (CAST)F_b_Vinch_K , -1, -1 }, // 1+3*2=7 //kj+++
{"hr_Vinch_K" , (CAST)F_hr_Vinch_K , -1, -1 }, // 1+3*1=4 //kj+++
{"Jr_Vinch_K" , (CAST)F_Jr_Vinch_K , -1, -1 }, // 1+3*1=4 //kj+++
{"h_Vinch_K" , (CAST)F_h_Vinch_K , -1, -1 }, // parameter is dimension {2},{3}, 3+3*2=9 //kj+++
{"dbdh_Vinch_K" , (CAST)F_dbdh_Vinch_K , -1, -1 }, // parameter is dimension {2},{3}, 1+3*2=7 //kj+++
{"dhdb_Vinch_K" , (CAST)F_dhdb_Vinch_K , -1, -1 }, // parameter is dimension {2},{3}, 1+3*2=7 //kj+++
// F_MultiHar
{"MHToTime" , (CAST)F_MHToTime , 0, 2 },
......
......@@ -69,6 +69,7 @@ void Cal_InitGalerkinTermOfFemEquation(struct EquationTerm * EquationTerm_P,
assDiag_done.clear();
//kj+++ // brutal approach: comment this to avoid auto-symmetrization of JacNL tensor with getdp (comment all this bloc)
// check for potentially symmetrical elementary matrix
FI->SymmetricalMatrix =
(EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ==
......@@ -94,6 +95,8 @@ void Cal_InitGalerkinTermOfFemEquation(struct EquationTerm * EquationTerm_P,
}
}
}
//*/ //kj+++
//FI->SymmetricalMatrix = 0; //kj+++ // brutal approach: uncomment this to avoid auto-symmetrization of JacNL tensor with getdp QQQ (uncomment this)
if (FI->SymmetricalMatrix) {
FI->Type_FormDof = FI->Type_FormEqu ;
......
......@@ -294,10 +294,17 @@ void F_dhdb_Vinch(F_ARG) ; // NOT USED FOR NOW (26/06/2016)
void F_dbdh_Vinch(F_ARG) ; // NOT USED FOR NOW (26/06/2016)
// kj+++
void F_Update_Jk (F_ARG) ; // NOT USED FOR NOW (26/06/2016)
void F_Update_Jk_sd (F_ARG) ; // NOT USED FOR NOW (26/06/2016)
void F_Update_Cell_K (F_ARG) ;
//Usefull Mathematical functions:
double Mul_VecVec_K(const double *v1, const double *v2);
void Mul_TensorVec_K(const double *M, const double *v, double *Mv, const int transpose_M);
void Mul_TensorSymTensorSym_K(double *A, double *B, double *C);
void Mul_TensorNonSymTensorNonSym_K(double *A, double *B, double *C);
void Mul_TensorNonSymTensorSym_K(double *A, double *B, double *C);
void Mul_TensorSymTensorNonSym_K(double *A, double *B, double *C);
void Inv_Tensor3x3_K (double *T, double *invT);
void Inv_TensorSym3x3_K (double *T, double *invT);
//Anhysteretic curve Characteristics:
double Lang(double nhr, double Ja, double ha);
double dLang(double nhr, double Ja, double ha);
double LangOverx(double nhr, double Ja, double ha);
......@@ -309,7 +316,6 @@ double dJanhy(double nhr,double Js, double alpha);
double Xanhy(double nhr,double Js, double alpha);
double dXanhy(double nhr,double Js, double alpha);
double IJanhy(double nhr, double Js, double alpha);
double Xanhy(double nhr,double Js, double alpha);
double InvJanhy(double nJ, double Js, double alpha);
double dInvJanhy(double nJ, double Js, double alpha) ;
......@@ -318,49 +324,91 @@ double dJanhy(double nJ, double Ja, double ha, double Jb, double hb) ;
double Xanhy(double nJ, double Ja, double ha, double Jb, double hb) ;
double dXanhy(double nJ, double Ja, double ha, double Jb, double hb) ;
double IJanhy(double nJ, double Ja, double ha, double Jb, double hb) ;
double Xanhy(double nJ, double Ja, double ha, double Jb, double hb) ;
double InvJanhy(double nJ, double Ja, double ha, double Jb, double hb) ;
double dInvJanhy_hr(double nhr, double Ja, double ha, double Jb, double hb);
double u_hr(double nhr, double Ja, double ha, double Jb, double hb) ;
double u_J(double nJ, double Js, double alpha);
double fct_omega(const double h[3], double Jk[3], const double Jkp[3], const double chi,
void Vector_Find_Jk_K (const double hrk[3],
const double Ja, const double ha, const double Jb, const double hb,
double Jk[3]);
void Vector_Find_hrk_K(const double Jk[3],
const double Ja, const double ha, const double Jb, const double hb,
double hrk[3]);
void Tensor_dJkdhrk_K(const double hr[3],
const double Ja, const double ha, const double Jb, const double hb,
double mutg[6]);
//Pseudo-Potential Functional Characteristics:
double fct_omega(const double h[3], const double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb);
void fct_d_omega (const double h[3], double Jk[3], const double Jkp[3], const double chi,
void fct_d_omega (const double h[3], const double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *d_omega);
void fct_dd_omega(const double h[3], const double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *ddfdJ2);
void fct_dd_omega_Num(const double h[3], const double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *ddfdJ2);
//void FullDiff_hi(double x, double chi, double ex[3]);
//double FullDiff_ff (double y, void *params);
//Usefull Functions for the Full Differential Approach:
void FullDiff_hi(double x, double chi, double ex[3]);
void FullDiff_xup(double x, double chi, double h[3], double xup[3]);
double FullDiff_ff_new (double y, void *params);
double FullDiff_ff (double y, void *params);
//Energy-Based Model - Vector Update:
void Vector_Update_Jk_K(const double h[3], double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb);
void Vector_Update_hr_K(const double h[3], double hr[3], const double hrp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb);
void Vector_Update_Simple_hr_K (const double h[3], double hr[3], const double hrp[3], const double chi) ;
void Vector_b_Vinch_K (const double h[3], double *Jk_all,
const double *Jkp_all, struct FunctionActive *D,double b[3] ) ;
void Vector_h_Vinch_K (const double b[3], double bc[3], double *Jk_all,
const double *Jkp_all, struct FunctionActive *D, double h[3] );
void Tensor_dbdh_Num(const double h[3], double *xk_all, const double *xkp_all,
struct FunctionActive *D,double *dbdh);
void Tensor_dhdb_Good_BFGS(const double dx[3],const double df[3],double *dhdb);
void Tensor_dbdh_Vinch_K ( const double h[3], double *Jk_all,
const double *Jkp_all, struct FunctionActive *D, double *dbdh);
void Tensor_dJkdh_Vinch_K(const double h[3], const double Jk[3], const double Jkp[3], const double chi,
//Energy-Based Model - Tensor Construction:
void Tensor_dJkdh_Var(const double h[3], const double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *dJkdh);
void Tensor_dJkdh_Var_Num(const double h[3], const double Jk[3], const double Jkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *dJkdh);
void Tensor_dJkdh_Diff_K( const double h[3], const double hrk[3], const double hrkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *dJkdh);
void Tensor_dhrkdh_Diff_ana(const double h[3], const double hrk[3], const double hrkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
const double dJkdhrk[6],
double *dhrkdh);
void Tensor_dhrkdh_Simple_ana(const double h[3], const double hrkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *dhrkdh);
void Tensor_dhrkdh_Num(const double h[3], const double xkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *dhrkdh);
void Tensor_dhrkdh_Simple_Num(const double h[3], const double xkp[3], const double chi,
double *dhrkdh); //NOT USE FOR NOW
void Tensor_dJkdh_Diff_K_old( const double h[3], const double hrk[3], const double hrkp[3], const double chi,
const double Ja, const double ha, const double Jb, const double hb,
double *dJkdh); //NOT USE FOR NOW
void Mul_TensorVec_K(const double *M, const double *v, double *Mv, const int transpose_M);
void Mul_TensorSymTensorSym_K(double *A, double *B, double *C);
double Mul_VecVec_K(const double *v1, const double *v2);
void Inv_Tensor3x3_K (double *T, double *invT);
void Inv_TensorSym3x3_K (double *T, double *invT);
void Tensor_dbdh_Vinch_K ( const double h[3], const double *Jk_all,
const double *Jkp_all, struct FunctionActive *D, double *dbdh);
void Tensor_dbdh_Num(const double h[3], double *xk_all, const double *xkp_all,
struct FunctionActive *D,double *dbdh);
void Tensor_dhdb_Good_BFGS(const double dx[3],const double df[3],double *dhdb);
//Energy-Based Model - GetDP Functions:
void F_Update_Cell_K (F_ARG) ;
void F_b_Vinch_K(F_ARG);
void F_hr_Vinch_K(F_ARG);
void F_Jr_Vinch_K(F_ARG);
void F_h_Vinch_K(F_ARG);
void F_dbdh_Vinch_K(F_ARG);
void F_dhdb_Vinch_K(F_ARG);
......
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -97,7 +97,7 @@ double CalcMaxErrorRatio(Resolution *Resolution_P,
else if (ErrorRatio > MaxErrorRatio)
MaxErrorRatio = ErrorRatio;
Current.Residual = ErrorRatio; //should be commented here QQQ?
//Current.Residual = ErrorRatio; //kj+++ (should be commented here, Current.ResidualN used instead and defined at the end of CalcMaxErrorRatio)
if (Message::GetVerbosity() > 5) {
Message::Info("IterativeLoopN: %s of %s error ratio from system %s: %.3g",
ILsystem.NormTypeString, ILsystem.NormOfString, DefineSystem_P->Name,
......@@ -143,7 +143,7 @@ double CalcMaxErrorRatio(Resolution *Resolution_P,
ErrorRatio);
}
}
Current.ResidualN = MaxErrorRatio ; //+++ Residual computed here for IterativeLoopN
Current.ResidualN = MaxErrorRatio ; //kj+++ Residual computed here for IterativeLoopN
return MaxErrorRatio;
}
......@@ -304,7 +304,8 @@ void Operation_IterativeLoopN(Resolution *Resolution_P,
if (Num_Iteration > NbrMaxIteration) {
Num_Iteration = NbrMaxIteration;
Flag_IterativeLoopConverged = 0;
Message::Info(3, "IterativeLoopN did NOT converge (%d iterations, error ratio %g)",
//Message::Info(3, "IterativeLoopN did NOT converge (%d iterations, error ratio %g)", //kj+++ (warning is better)
Message::Warning("IterativeLoopN did NOT converge (%d iterations, error ratio %g)",
(int)Current.Iteration, MaxErrorRatio);
}
Current.Iteration = Save_Iteration ;
......
......@@ -1271,7 +1271,7 @@ void Treatment_Operation(struct Resolution * Resolution_P,
LinAlg_AssembleVector(&DofData_P->b) ;
}
//LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &DofData_P->dx); //In prevision to build 'dx' in the following (needed for "IterativeLoopPro") QQQ?
LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &DofData_P->dx); //kj+++ In prevision to build 'dx' in the following (needed for "IterativeLoopPro")
if(!again){
LinAlg_Solve(&DofData_P->A, &DofData_P->b, &DofData_P->Solver,
......@@ -1286,7 +1286,7 @@ void Treatment_Operation(struct Resolution * Resolution_P,
(Operation_P->Flag < 0) ? 0 : Operation_P->Flag) ;
}
//LinAlg_SubVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx, &DofData_P->dx) ; //In order to build 'dx' (needed for "IterativeLoopPro") QQQ?
LinAlg_SubVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx, &DofData_P->dx) ; //kj+++ In order to build 'dx' (needed for "IterativeLoopPro")
#ifdef TIMER
double timer = MPI_Wtime() - tstart;
printf("Proc %d, time spent in %s %.16g\n", again ? "SolveAgain" : "Solve",
......@@ -1375,18 +1375,19 @@ void Treatment_Operation(struct Resolution * Resolution_P,
else
LinAlg_SolveAgain(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver, &DofData_P->dx) ;
//............... The following lines (*) are not needed here because ..............
// Current.Residual is only needed when Flag_IterativeLoopN==0 QQQ?
Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0, &MeanError) ;
// NormType: 1=LINFNORM, 2=L1NORM, 3=MEANL1NORM, 4=L2NORM, 5=MEANL2NORM
// Closest behaviour to old function Cal_SolutionError with MEANL2NORM
// Cal_SolutionErrorRatio(&DofData_P->dx, &DofData_P->CurrentSolution->x,
// reltol, abstol, MEANL2NORM, &MeanError) ;
//LinAlg_VectorNorm2(&DofData_P->dx, &MeanError);
Current.Residual = MeanError; // NB: Residual computed for classical IterativeLoop using SolveJac
//.................................................................................
if(!Flag_IterativeLoopN){
//............... previous lines (*) should be placed here QQQ? ..................
//kj+++ (the following lines were outside the condition (!Flag_IterativeLoopN) but
// Current.Residual is only needed when Flag_IterativeLoopN==0)
//.........................................................................................
Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0, &MeanError) ;
// NormType: 1=LINFNORM, 2=L1NORM, 3=MEANL1NORM, 4=L2NORM, 5=MEANL2NORM
// Closest behaviour to old function Cal_SolutionError with MEANL2NORM
// Cal_SolutionErrorRatio(&DofData_P->dx, &DofData_P->CurrentSolution->x,
// reltol, abstol, MEANL2NORM, &MeanError) ;
//LinAlg_VectorNorm2(&DofData_P->dx, &MeanError);
Current.Residual = MeanError; // NB: Residual computed for classical IterativeLoop using SolveJac
//.........................................................................................
if (MeanError != MeanError){
Message::Warning("No valid solution found (NaN or Inf)!");
}
......@@ -1457,7 +1458,7 @@ void Treatment_Operation(struct Resolution * Resolution_P,
/* MHBilinear-terms do not contribute to the RHS and residual, and are thus disregarded */
Error_Prev = 1e99 ; Frelax_Opt = 1. ;
//if(Current.Iteration==1) Current.Residual_Iter1=1.0; //to compute a relative residual (relative to residual at iter 1) in SolveJacAdapt //QQQ?
//if(Current.Iteration==1) Current.Residual_Iter1=1.0; //kj+++ to compute a relative residual (relative to residual at iter 1) in SolveJacAdapt (not necessary)
if (!(NbrSteps_relax = List_Nbr(Operation_P->Case.SolveJac_AdaptRelax.Factor_L))){
Message::Error("No factors provided for Adaptive Relaxation");
......@@ -1491,7 +1492,7 @@ void Treatment_Operation(struct Resolution * Resolution_P,
LinAlg_VectorNorm2(&DofData_P->res, &Norm);
LinAlg_GetVectorSize(&DofData_P->res, &N);
Norm /= (double)N;
//Norm /= Current.Residual_Iter1; //to compute a relative residual (relative to residual at iter 1) in SolveJacAdapt //QQQ?
//Norm /= Current.Residual_Iter1; //kj+++ to compute a relative residual (relative to residual at iter 1) in SolveJacAdapt (not necessary)
Current.Residual = Norm;
if(Message::GetVerbosity() == 10)
......@@ -1514,11 +1515,11 @@ void Treatment_Operation(struct Resolution * Resolution_P,
Frelax_Opt, &DofData_P->CurrentSolution->x);
MeanError = Error_Prev ;
Current.RelaxFac = Frelax_Opt; // +++
//Current.Residual = MeanError; // QQQ? Residual computed here with SolveJacAdapt (usefull to test stop criterion in classical IterativeLoop then)
Current.RelaxFac = Frelax_Opt; //kj+++
Current.Residual = MeanError; //kj+++ Residual computed here with SolveJacAdapt (usefull to test stop criterion in classical IterativeLoop then)
Message::Info("%3ld Nonlinear Residual norm %14.12e (optimal relaxation factor = %f)",
(int)Current.Iteration, MeanError, Frelax_Opt);
//if(Current.Iteration==1) Current.Residual_Iter1=MeanError; //to compute a relative residual (relative to residual at iter 1) in SolveJacAdapt //QQQ?
//if(Current.Iteration==1) Current.Residual_Iter1=MeanError; //kj+++ to compute a relative residual (relative to residual at iter 1) in SolveJacAdapt (not necessary)
if(Message::GetProgressMeterStep() > 0 && Message::GetProgressMeterStep() < 100)
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Residual",
std::vector<double>(1, MeanError));
......@@ -2596,7 +2597,7 @@ void Treatment_Operation(struct Resolution * Resolution_P,
Save_Iteration = Current.Iteration ;
//Current.Residual_Iter1=1.0; //to compute a relative residual (relative to residual at iter 1) QQQ?
//Current.Residual_Iter1=1.0; //kj+++ to compute a relative residual (relative to residual at iter 1) (not necessary)
// abstol = Operation_P->Case.IterativeLoop.Criterion ;
// reltol = Operation_P->Case.IterativeLoop.Criterion/100 ;
......@@ -2624,13 +2625,13 @@ void Treatment_Operation(struct Resolution * Resolution_P,
// TransferSolution of a nonlinear resolution
Treatment_Operation(Resolution_P, Operation_P->Case.IterativeLoop.Operation,
DofData_P0, GeoData_P0, Resolution2_P, DofData2_P0) ;
//if (Current.Iteration==1) Current.Residual_Iter1=Current.RelativeDifference; //to compute a relative residual (relative to residual at iter 1) QQQ?
//if (Current.Iteration==1) Current.Residual_Iter1=Current.RelativeDifference; //kj+++ to compute a relative residual (relative to residual at iter 1) (not necessary)
//NB: Current.RelativeDifference is what is used for classical IterativeLoop stop criterion
//NB: In SolveJac: (Current.RelativeDifference+=Current.Residual)
//NB: In SolveJacAdapt: (Current.RelativeDifference=Current.Residual)
if ( (Current.RelativeDifference <= Operation_P->Case.IterativeLoop.Criterion) ||
//(Current.RelativeDifference/Current.Residual_Iter1 <= Operation_P->Case.IterativeLoop.Criterion) || //to compute a relative residual (relative to residual at iter 1) QQQ?
//(Current.RelativeDifference/Current.Residual_Iter1 <= Operation_P->Case.IterativeLoop.Criterion) || //kj+++ to compute a relative residual (relative to residual at iter 1) (not necessary)
(Current.RelativeDifference != Current.RelativeDifference) ) // NaN or Inf
break ;
......@@ -2655,16 +2656,16 @@ void Treatment_Operation(struct Resolution * Resolution_P,
Num_Iteration, Num_Iteration > 1 ? "s" : "",
Current.RelativeDifference);
}
/* QQQ?
/* kj+++ to compute a relative residual (relative to residual at iter 1) (not necessary)
Message::Info(3, "IterativeLoop did NOT converge (%d iterations, residual %g)\n rel %g",
Num_Iteration, Current.RelativeDifference, Current.RelativeDifference/Current.Residual_Iter1_kj);
Num_Iteration, Current.RelativeDifference, Current.RelativeDifference/Current.Residual_Iter1);
// Either it has reached the max num of iterations or a NaN at a given iteration
Num_Iteration = Operation_P->Case.IterativeLoop.NbrMaxIteration ;
}
else{
Message::Info(3, "IterativeLoop converged (%d iteration%s, residual %g)\n rel %g",
Num_Iteration, Num_Iteration > 1 ? "s" : "",
Current.RelativeDifference, Current.RelativeDifference/Current.Residual_Iter1_kj);
Current.RelativeDifference, Current.RelativeDifference/Current.Residual_Iter1);
}
*/
Current.Iteration = Save_Iteration ;
......
......@@ -4,15 +4,28 @@ Group{
Function{
GetItLoopInfo[]=Tensor[ $iter,
$res,
$res/$res0,
$inc,
$relax,
$relaxcount,
$relaxcounttot,
$syscount,
$dnccount];
If(Flag_NLRes==NLRES_ITERATIVELOOPPRO)
GetItLoopInfo[]=Tensor[ $iter,
$res,
$res/$res0,
$inc,
$relax,
$relaxcount,
$relaxcounttot,
$syscount,
$dnccount];
Else
GetItLoopInfo[]=Tensor[ $iter,
$res,
$resL,
$resN,
$relax,
$relaxcount,
$relaxcounttot,
$syscount,
$dnccount];
EndIf
//--------------------------------------
// M Y . I T L O O P P R O -------------
......@@ -63,9 +76,8 @@ Macro MyItLoopPro
GetNormIncrement[A, $inc]; //Check the new increment
GetNormSolution[A, $sol]; //Check the new solution
Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, $relaxcount}, Format
"*** ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relaxopt = %1.3g, relaxcount=%02g ***"];
PostOperation[PostOpItLoop_a];
Call DisplayRunTimeVar;
}
Return
//...............................................................
......@@ -194,6 +206,26 @@ Macro MySolveJac_AdaptRelax
$res=$res_Prev];
Return
//...............................................................
//--------------------------------------------------------
// D I S P L A Y . R U N . T I M E . V A R ---------------
//--------------------------------------------------------
Macro DisplayRunTimeVar
If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
Print[{$TimeStep, $iter, $res, $resL, $resN, $relax, $relaxcount}, Format
"*** ts=%05g, it=%03g, res_used=%1.5e, (res = %1.5e, resN=%1.5e), relaxopt = %1.3g, relaxcount=%02g ***"];
EndIf
If(Flag_NLRes==NLRES_ITERATIVELOOPPRO)
Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, $relaxcount}, Format
"*** ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relaxopt = %1.3g, relaxcount=%02g ***"];
EndIf
PostOperation[PostOpItLoop_a];
Return
//...............................................................
}
// To Print Info about Convergence
......
......@@ -192,12 +192,12 @@ Integration {
{
Type Gauss ;
Case {
{ GeoElement Triangle ; NumberOfPoints 1 ; }
{ GeoElement Triangle ; NumberOfPoints 1 ; }
{ GeoElement Quadrangle ; NumberOfPoints 1 ; }
{ GeoElement Tetrahedron ; NumberOfPoints 1 ; }
{ GeoElement Hexahedron ; NumberOfPoints 1 ; }
{ GeoElement Prism ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 1 ; }
{ GeoElement Line ; NumberOfPoints 1 ; }
}
}
}
......@@ -577,9 +577,9 @@ Formulation {
//-----------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------------
If(Flag_NLRes==NLRES_ITERATIVELOOPPRO)
Include "MyItLoopPro.pro";
EndIf
//If(Flag_NLRes==NLRES_ITERATIVELOOPPRO)
Include "MyItLoopPro.pro";
//EndIf
//-----------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------------
......@@ -632,6 +632,7 @@ Resolution {
Evaluate[$relaxcount=0];
Evaluate[$relaxcounttot=0];
Evaluate[$dnccount=0];
Evaluate[$itertot=0];
//Save TimeStep 0
SaveSolution[A];
......@@ -644,6 +645,7 @@ Resolution {
EndIf
If(Flag_NL) // NonLinear Case ...
Evaluate[$relax=relax_max];
If(Flag_NLRes==NLRES_ITERATIVELOOP)
// I T E R A T I V E . L O O P ...................................
IterativeLoop[ Nb_max_iter, stop_criterion, RF_tuned[]] {
......@@ -654,23 +656,45 @@ Resolution {
Else
SolveJac_AdaptRelax[A, List[RelaxFac_Lin],TestAllFactors ] ;
Evaluate[$relaxcount=$NbrTestedFac];
EndIf
EndIf
Evaluate[$res = $Residual, $resL = $Residual,$resN = $ResidualN, $iter = $Iteration];
Evaluate[$syscount = $syscount + 1 ];
Evaluate[$relaxcounttot=$relaxcounttot+$relaxcount];
Evaluate[$relax= $RelaxFac];
Call DisplayRunTimeVar;
}
Evaluate[$syscount = $syscount + 1 ];
Evaluate[$relaxcounttot=$relaxcounttot+$relaxcount];
Evaluate[$relax= $RelaxFac];
//...............................................................
EndIf
If(Flag_NLRes==NLRES_ITERATIVELOOPN)
// I T E R A T I V E . L O O P . N ..............................
//:::::::::::::::::::::::::::::::::
// INIT for h_only case::::::::::::
/*
GenerateJac[A] ;
If (!Flag_AdaptRelax)
SolveJac[A] ;
Else
SolveJac_AdaptRelax[A, List[RelaxFac_Lin],TestAllFactors ] ;
EndIf
//*/
//:::::::::::::::::::::::::::::::::
IterativeLoopN[ Nb_max_iter, RF_tuned[],
//*****Choose between one of the 3 following possibilities:*****
//System { { A , Reltol_Mag, Abstol_Mag, Residual MeanL2Norm }} ]{ //1)
//PostOperation { { a_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //2)
//PostOperation { { b_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //3)
PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //4)
PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //4) //Need the above "INIT" for square with EnergHyst or JA because h_only = 0 at iter 1
//**************************************************************
Evaluate[$res = $ResidualN, $resL = $Residual,$resN = $ResidualN,$iter = $Iteration-1];
Test[ $iter>0]{
Call DisplayRunTimeVar;
}
GenerateJac[A] ;
If (!Flag_AdaptRelax)
SolveJac[A] ;
......@@ -679,10 +703,14 @@ Resolution {
SolveJac_AdaptRelax[A, List[RelaxFac_Lin],TestAllFactors ] ;
Evaluate[$relaxcount=$NbrTestedFac];
EndIf
Evaluate[$syscount = $syscount + 1 ];
Evaluate[$relaxcounttot=$relaxcounttot+$relaxcount];
Evaluate[$relax= $RelaxFac];
}
Evaluate[$syscount = $syscount + 1 ];
Evaluate[$relaxcounttot=$relaxcounttot+$relaxcount];
Evaluate[$relax= $RelaxFac];
Evaluate[$res = $ResidualN, $resL = $Residual, $resN = $ResidualN, $iter = $iter+1];
Call DisplayRunTimeVar;
//...............................................................
EndIf
......@@ -691,14 +719,35 @@ Resolution {
Call MyItLoopPro;
//...............................................................
EndIf
EndIf
Evaluate[$itertot = $itertot+$iter];
Test[ (Flag_NLRes==NLRES_ITERATIVELOOPPRO && $iter == Nb_max_iter) ||
(Flag_NLRes==NLRES_ITERATIVELOOP && $res > Abstol_Mag ) ||
(Flag_NLRes==NLRES_ITERATIVELOOPN && $res > 1. ) ]{
Print[{$Time, $TimeStep, $DTime, $relax, $iter}, Format
"*** Non convergence at Theta Time = %g s (TimeStep %g, DTime %g, relax %g) (iter %g)"];
Print["--> Save bad solution and move to the next time step"];
Evaluate[$dnccount=$dnccount+1];
}
EndIf
SaveSolution[A];
PostOperation[Get_LocalFields] ;
If (Flag_LiveLocalPostOp)
PostOperation[Get_LocalFields] ;
EndIf
//Test[ $TimeStep > 1 ]{
PostOperation[Get_GlobalQuantities];
If (Flag_LiveGlobalPostOp)
PostOperation[Get_GlobalQuantities];
EndIf
//}
}
Print["--------------------------------------------------------------"];
Print[{$syscount}, Format "Total number of linear systems solved: %g"];
Print[{$relaxcounttot}, Format "Total number of relaxation factor tested: %g"];
Print[{$itertot/$TimeStep}, Format "Mean number of NR iter at FE level: %g"];
Print[{$dnccount}, Format "Total number of non converged TimeStep: %g"];
Print["--------------------------------------------------------------"];
EndIf // End Flag_AnalysisType==AN_TIME (Time domain)
}
}
......
dim00 = (Flag_3Dmodel==1)?3:2;
N00 = 3;
FLAG_TANORLANG00 = 1;
FLAG_VARORDIFF00 = 2;
FLAG_VARORDIFF00 = 1;
FLAG_INVMETHOD00 = 1;
FLAG_ROOTFINDING1D00 = 0;
FLAG_MINMETHOD00 = 1;
......@@ -14,4 +14,5 @@ TOLERANCE_OM00 = 1.e-11;
MAX_ITER_OM00 = 700;
FLAG_ANA00 = 1;
TOLERANCE_NJ00 = 1.e-3;
DELTA_000 = 1.e-0;
\ No newline at end of file
DELTA_000 = 1.e-0;
FLAG_HOMO00 = 0;
\ No newline at end of file
......@@ -50,6 +50,16 @@ FLAG_MINMETHOD = 6 --> steepest descent (gsl)
*/
FLAG_WARNING = FLAG_WARNING00 ; // SHOW WARNING
/*
#define FLAG_WARNING_INFO_INV 1
#define FLAG_WARNING_INFO_ROOTFINDING 2
#define FLAG_WARNING_INFO_MIN 3
#define FLAG_WARNING_DISP_INV 11
#define FLAG_WARNING_DISP_ROOTFINDING 12
#define FLAG_WARNING_STOP_INV 101
#define FLAG_WARNING_STOP_ROOTFINDING 102
#define FLAG_WARNING_ITER 100
*/
TOLERANCE_JS = TOLERANCE_JS00; // SENSITIVE_PARAM (1.e-3) // 1.e-4
TOLERANCE_0 = TOLERANCE_000; // SENSITIVE_PARAM (1.e-7)
TOLERANCE_NR = TOLERANCE_NR00; // SENSITIVE_PARAM (1.e-7) // 1.e-8 needed for diff with NR,1.e-5
......@@ -61,22 +71,51 @@ FLAG_MINMETHOD = 6 --> steepest descent (gsl)
// 1.e-3 for VinchT.pro & transfo.pro)
DELTA_0 = DELTA_000; // SENSITIVE_PARAM (1.e-3 for square;
// 1.e0 for VinchT & transfo)
FLAG_HOMO = FLAG_HOMO00;
If(FLAG_TANORLANG==1)
// Material M25050A
Ja = 1.22;
ha = 65;
Jb = 0.;
hb = 0.;
EndIf
If(FLAG_TANORLANG==2)
/*
// Material M23535A
Ja = 0.595;
ha = 4100;
Jb = 1.375;
hb = 17.5;
/*/
/* //team32_RD_z_anhyshift (abstract ISEM 2017)
Ja = 1.03838;
ha = 16.9403;
Jb = 0.600978;
hb = 240.578;
//*/
/* //Double Langevin for team32_RD_z_anhyFab (new try)
Ja=1.03643 ;
ha=14.0374 ;
Jb=0.549518 ;
hb=188.516;
//*/
///* //FH(xini=400) // NEW 1/9/2017
Ja = 0.792234;
ha = 9.082095;
Jb = 0.790993;
hb = 137.121351;
//*/
/* //TD
Ja = 1.08427393613
ha = 29.1025536684
Jb = 0.572282990517
hb = 399.474553851
*/
EndIf
If(N==3)
///*
// Material M25050A
Js_1 = 0.11;
Js_2 = 0.8;
Js_3 = 0.31;
......@@ -86,6 +125,46 @@ FLAG_MINMETHOD = 6 --> steepest descent (gsl)
w_1 = 0.090163934426230;
w_2 = 0.655737704918033;
w_3 = 0.254098360655738;
/*/
/* //team32_RD_z_anhyshift (abstract ISEM 2017)
chi_1 = 0*11.6413;
chi_2 = 50.1413;
chi_3 = 107.518;
w_1 = 0.173524;
w_2 = 0.591535;
w_3 = 0.234941;
Js_1 = w_1*(Ja+Jb);
Js_2 = w_2*(Ja+Jb);
Js_3 = w_3*(Ja+Jb);
//*/
/* //FH(xini=400) // THIS IS PROBLEMATIC (CONVERGENCE PROBLEM ok now) (new try)
chi_1 = 0.;
chi_2 = 53.779614004;
chi_3 = 233.384447699;
w_1 = 0.1048; //+0.1 to converge
w_2 = 0.817943263499; // -0.1 to converge
w_3 = 0.0772371911006;
Js_1 = w_1*(Ja+Jb);
Js_2 = w_2*(Ja+Jb);
Js_3 = w_3*(Ja+Jb);
//*/
/* //FH(xini=400) // NEW 1/9/2017 (chamonix)
chi_1 = 0.;
chi_2 = 53.789872;
chi_3 = 233.910892;
w_1 = 0.0280702473932;
w_2 = 0.896565361832;
w_3 = 0.0753643907746;
Js_1 = w_1*(Ja+Jb);
Js_2 = w_2*(Ja+Jb);
Js_3 = w_3*(Ja+Jb);
//*/
/*//FH(xini=400) TD // NEW 1/9/2017 (chamonix)
0 0.0565084014449 0.0
1 0.839372763961 65.418004
2 0.104118834594 246.564466
*/
EndIf
If(N==1)
Js_1 = 1;
......@@ -98,12 +177,64 @@ FLAG_MINMETHOD = 6 --> steepest descent (gsl)
w_2 = 0.;
w_3 = 0.;
EndIf
If(N==5) //FH(xini=400) // NEW 1/9/2017
chi_1 = 0.;