Commit fcecb101 authored by Kevin Jacques's avatar Kevin Jacques

fix F.h + kj modif on Residual and dx applied in SolvingOperations.cpp and…

fix F.h + kj modif on Residual and dx applied in SolvingOperations.cpp and Operation_IterativeLoopN.cpp, kj modif for avoiding auto symmetrisation only suggested in Cal_GalerkinTermOfFemEquation.cpp
parent e9fb8e74
Pipeline #1126 passed with stage
in 44 seconds
......@@ -69,7 +69,7 @@ void Cal_InitGalerkinTermOfFemEquation(struct EquationTerm * EquationTerm_P,
assDiag_done.clear();
//QQQ??? // brutal approach: comment this to avoid auto-symmetrization of JacNL tensor with getdp QQQ (comment all this)
//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 ==
......@@ -95,8 +95,8 @@ void Cal_InitGalerkinTermOfFemEquation(struct EquationTerm * EquationTerm_P,
}
}
}
//*/ //QQQ???
//FI->SymmetricalMatrix = 0; //QQQ??? // brutal approach: uncomment this to avoid auto-symmetrization of JacNL tensor with getdp QQQ (uncomment this)
//*/ //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);
......
......@@ -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,12 +304,9 @@ 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)",
(int)Current.Iteration, MaxErrorRatio);
/* //QQQ??? (uncomment the following and comment the previous Message)
//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 ;
Flag_IterativeLoopN = 0;
......
......@@ -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??? (uncomment this)
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? (uncomment this)
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 (move this)
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 (here)..................
//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? (not necessary)
//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? (not necessary)
//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) //QQQ (uncomment this)
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? (not necessary)
//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? (not necessary)
//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? (not necessary)
//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? (not necessary)
//(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? (not necessary)
/* 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 ;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment