"Fossies" - the Fresh Open Source Software Archive

Member "getdp-3.1.0-source/Kernel/Operation_IterativeTimeReduction.cpp" (29 Dec 2018, 19715 Bytes) of package /linux/privat/getdp-3.1.0-source.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Operation_IterativeTimeReduction.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 3.0.4_vs_3.1.0.

    1 // GetDP - Copyright (C) 1997-2019 P. Dular and C. Geuzaine, University of Liege
    2 //
    3 // See the LICENSE.txt file for license information. Please report all
    4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
    5 
    6 #include <stdio.h>
    7 #include <stdlib.h>
    8 #include <math.h>
    9 #include "ProData.h"
   10 #include "DofData.h"
   11 #include "SolvingOperations.h"
   12 #include "Message.h"
   13 #include "MallocUtils.h"
   14 #include "Cal_Quantity.h"
   15 #include "Get_DofOfElement.h"
   16 
   17 #define TWO_PI     6.2831853071795865
   18 
   19 extern struct Problem Problem_S ;
   20 extern struct CurrentData Current ;
   21 
   22 extern int  Flag_NextThetaFixed ;
   23 
   24 /* ------------------------------------------------------------------------ */
   25 /*  C a l _ S o l u t i o n E r r o r X                                     */
   26 /* ------------------------------------------------------------------------ */
   27 
   28 void Cal_SolutionErrorX(int Nbr, double * xNew, double * x, double * MeanError)
   29 {
   30   int i;
   31   double errsqr = 0., xmoy = 0., dxmoy = 0., tol ;
   32 
   33   if(0 && gSCALAR_SIZE == 2)
   34     Message::Error("FIXME: Cal_SolutionErrorX might return strange results"
   35                    " in complex arithmetic");
   36 
   37   for (i = 0 ; i < Nbr ; i++) {
   38     xmoy  += fabs( x[i])/(double)Nbr ;
   39     dxmoy += fabs(xNew[i]-x[i])/(double)Nbr ;
   40   }
   41 
   42   if (xmoy > 1.e-30) {
   43     tol = xmoy*1.e-10 ;
   44 
   45     for (i = 0 ; i < Nbr ; i++)
   46       if ( fabs(x[i]) > tol ) errsqr += fabs((xNew[i]-x[i]) / x[i]) ;
   47       else                    errsqr += fabs(xNew[i]-x[i]) ;
   48 
   49     *MeanError = errsqr / (double)Nbr ;
   50   }
   51   else
   52     if (dxmoy > 1.e-30)  *MeanError = 1. ;
   53     else                 *MeanError = 0. ;
   54 }
   55 
   56 /* ------------------------------------------------------------------------ */
   57 /*  C a l _ C o m p a r e G l o b a l Q u a n t i t y                       */
   58 /* ------------------------------------------------------------------------ */
   59 
   60 #define COMPARE_CHANGE        1
   61 #define COMPARE_CONVERGENCE   2
   62 
   63 void  Cal_CompareGlobalQuantity(struct Operation * Operation_P,
   64                 int Type_Analyse, int * Type_ChangeOfState,
   65                 int * FlagIndex, int Flag_First)
   66 {
   67   List_T  *Region_L = NULL ;
   68   int      i, Nbr_Region = 0, Num_Region ;
   69   int      Nbr_ChangeOfState, i_COS ;
   70 
   71   struct ChangeOfState     *ChangeOfState_P = NULL;
   72   struct Formulation       *Formulation_P ;
   73   struct FunctionSpace     *FunctionSpace_P ;
   74   struct GlobalQuantity    *GlobalQuantity_P ;
   75   struct DefineQuantity    *DefineQuantity_P ;
   76   struct QuantityStorage   QuantityStorage_S ;
   77 
   78   double  Val0_Dof, Val1_Dof ;
   79   double  *val0 = NULL, *val1 = NULL, MeanError, v0, v1 ;
   80   struct  Value  Value ;
   81 
   82   double  Val1_E, Val0_E, Val_S, Val0_Ref, Val1_Ref, v_fz, v_k, v_ke, v_sat ;
   83   double  Save_Time ;
   84 
   85   if(0 && gSCALAR_SIZE == 2)
   86     Message::Error("FIXME: Cal_CompareGlobalQuantity might return strange results"
   87                    " in complex arithmetic");
   88 
   89   /* test */
   90   v_k  = 1./27.2836 ;  v_ke = 18.518519 ;
   91   v_fz = 1. / (5.e3 * 1.5e-9) ;  v_sat = 6. ;
   92 
   93 
   94   *Type_ChangeOfState = CHANGEOFSTATE_NOCHANGE ;
   95 
   96   Nbr_ChangeOfState =
   97     List_Nbr(Operation_P->Case.IterativeTimeReduction.ChangeOfState) ;
   98 
   99   for (i_COS = 0 ; i_COS < Nbr_ChangeOfState ; i_COS++) {
  100 
  101     ChangeOfState_P = (struct ChangeOfState *)
  102       List_Pointer(Operation_P->Case.IterativeTimeReduction.ChangeOfState, i_COS) ;
  103 
  104     Region_L =
  105       ((struct Group *)List_Pointer(Problem_S.Group, ChangeOfState_P->InIndex))
  106       ->InitialList ;
  107     List_Sort(Region_L, fcmp_int) ;
  108 
  109     Nbr_Region = List_Nbr(Region_L) ;
  110 
  111     if (Nbr_Region > 0) {
  112 
  113       Formulation_P = (struct Formulation *)
  114     List_Pointer(Problem_S.Formulation, ChangeOfState_P->FormulationIndex) ;
  115       DefineQuantity_P = (struct DefineQuantity *)
  116     List_Pointer(Formulation_P->DefineQuantity, ChangeOfState_P->QuantityIndex) ;
  117 
  118       QuantityStorage_S.FunctionSpace = FunctionSpace_P =
  119     (struct FunctionSpace*)List_Pointer(Problem_S.FunctionSpace,
  120                         DefineQuantity_P->FunctionSpaceIndex) ;
  121       GlobalQuantity_P = (struct GlobalQuantity*)
  122     List_Pointer(FunctionSpace_P->GlobalQuantity,
  123              *(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
  124 
  125       if (!ChangeOfState_P->ActiveList[0]) {
  126     ChangeOfState_P->ActiveList[0] =(double *)Malloc(Nbr_Region * sizeof(double)) ;
  127     ChangeOfState_P->ActiveList[1] =(double *)Malloc(Nbr_Region * sizeof(double)) ;
  128       }
  129       val0 = ChangeOfState_P->ActiveList[0] ; val1 = ChangeOfState_P->ActiveList[1] ;
  130 
  131       /* debug */
  132 
  133       if (Type_Analyse == 999 && i_COS == 0 &&
  134       ChangeOfState_P->Type == CHANGEOFSTATE_CHANGEREFERENCE2) {
  135 
  136     List_Read(Region_L, 0, &Num_Region) ;
  137     Current.Region = Num_Region ;
  138 
  139     Get_DofOfRegion(Current.Region, GlobalQuantity_P,
  140             FunctionSpace_P, &QuantityStorage_S) ;
  141 
  142 
  143     QuantityStorage_S.FunctionSpace->DofData->CurrentSolution -- ;
  144     Save_Time = Current.Time ;
  145     Current.Time =
  146       QuantityStorage_S.FunctionSpace->DofData->CurrentSolution->Time ;
  147     Dof_GetRealDofValue(QuantityStorage_S.FunctionSpace->DofData,
  148                 QuantityStorage_S.BasisFunction[0].Dof, &Val0_Dof) ;
  149 
  150     Get_ValueOfExpressionByIndex(ChangeOfState_P->ExpressionIndex,
  151                      NULL, 0., 0., 0., &Value) ;
  152     Val0_Ref = Value.Val[0] ;
  153 
  154     Current.Time = Save_Time ;
  155     QuantityStorage_S.FunctionSpace->DofData->CurrentSolution ++ ;
  156 
  157     Dof_GetRealDofValue(QuantityStorage_S.FunctionSpace->DofData,
  158                 QuantityStorage_S.BasisFunction[0].Dof, &Val1_Dof) ;
  159 
  160     Get_ValueOfExpressionByIndex(ChangeOfState_P->ExpressionIndex,
  161                      NULL, 0., 0., 0., &Value) ;
  162     Val1_Ref = Value.Val[0] ;
  163 
  164     Val1_E = (Val1_Ref - v_k * Val1_Dof) * v_ke ;
  165     Val0_E = (Val0_Ref - v_k * Val0_Dof) * v_ke ;
  166 
  167     Val_S = Val1_E + (Val1_E-Val0_E)/Current.DTime / (TWO_PI*v_fz) ;
  168     /*
  169       fprintf(FilePWM, "%.16g %g %g", Current.Time, Val1_E, Val_S) ;
  170     */
  171     Val_S += Val1_Ref ;
  172     if      (Val_S >  v_sat)  Val_S =  v_sat ;
  173     else if (Val_S < -v_sat)  Val_S = -v_sat ;
  174     /*
  175       fprintf(FilePWM, " %g %g\n", Val_S,
  176       ((struct Expression *)
  177       List_Pointer(Problem_S.Expression, ChangeOfState_P->FlagIndex))
  178       ->Case.Constant
  179       ) ;
  180 
  181       fflush(FilePWM) ;
  182     */
  183     break ;
  184       }
  185       /* else if (Type_Analyse == 999 && i_COS > 0) break ; */
  186 
  187       /* ----- */
  188 
  189       /*  C a l c u l   v a l e u r s   . . .  */
  190 
  191       for (i = 0 ; i < Nbr_Region ; i++) {
  192 
  193     List_Read(Region_L, i, &Num_Region) ;
  194     Current.Region = Num_Region ;
  195 
  196     if (DefineQuantity_P->Type == GLOBALQUANTITY) {
  197       Get_DofOfRegion(Current.Region, GlobalQuantity_P,
  198               FunctionSpace_P, &QuantityStorage_S) ;
  199       switch (Type_Analyse) {
  200       case COMPARE_CHANGE :  /* Compare values at times t-dt and t */
  201         Dof_GetRealDofValue(QuantityStorage_S.FunctionSpace->DofData,
  202                 QuantityStorage_S.BasisFunction[0].Dof, &Val1_Dof) ;
  203 
  204         switch (ChangeOfState_P->Type) {
  205         case CHANGEOFSTATE_CHANGESIGN :
  206         case CHANGEOFSTATE_CHANGELEVEL :
  207           QuantityStorage_S.FunctionSpace->DofData->CurrentSolution -- ;
  208           Dof_GetRealDofValue(QuantityStorage_S.FunctionSpace->DofData,
  209                   QuantityStorage_S.BasisFunction[0].Dof, &Val0_Dof) ;
  210           QuantityStorage_S.FunctionSpace->DofData->CurrentSolution ++ ;
  211           break ;
  212 
  213         case CHANGEOFSTATE_CHANGEREFERENCE :
  214           Get_ValueOfExpressionByIndex(ChangeOfState_P->ExpressionIndex,
  215                        NULL, 0., 0., 0., &Value) ;
  216           Val0_Dof = Value.Val[0] ;
  217 
  218           break ;
  219 
  220         case CHANGEOFSTATE_CHANGEREFERENCE2 :
  221           QuantityStorage_S.FunctionSpace->DofData->CurrentSolution -- ;
  222           Save_Time = Current.Time ;
  223           Current.Time =
  224         QuantityStorage_S.FunctionSpace->DofData->CurrentSolution->Time ;
  225           Dof_GetRealDofValue(QuantityStorage_S.FunctionSpace->DofData,
  226                   QuantityStorage_S.BasisFunction[0].Dof, &Val0_Dof) ;
  227           Get_ValueOfExpressionByIndex(ChangeOfState_P->ExpressionIndex,
  228                        NULL, 0., 0., 0., &Value) ;
  229           Val0_Ref = Value.Val[0] ;
  230           Current.Time = Save_Time ;
  231           QuantityStorage_S.FunctionSpace->DofData->CurrentSolution ++ ;
  232 
  233           Get_ValueOfExpressionByIndex(ChangeOfState_P->ExpressionIndex,
  234                        NULL, 0., 0., 0., &Value) ;
  235           Val1_Ref = Value.Val[0] ;
  236 
  237           Val1_E = (Val1_Ref - v_k * Val1_Dof) * v_ke ;
  238           Val0_E = (Val0_Ref - v_k * Val0_Dof) * v_ke ;
  239 
  240           Val_S = Val1_E + (Val1_E-Val0_E)/Current.DTime / (TWO_PI*v_fz) ;
  241           Val_S += Val1_Ref ;
  242           if      (Val_S >  v_sat)  Val_S =  v_sat ;
  243           else if (Val_S < -v_sat)  Val_S = -v_sat ;
  244 
  245           Val1_Dof = Val_S ;
  246 
  247           Get_ValueOfExpressionByIndex(ChangeOfState_P->ExpressionIndex2,
  248                        NULL, 0., 0., 0., &Value) ;
  249           Val0_Dof = Value.Val[0] ;
  250 
  251           break ;
  252         }
  253         break ;
  254 
  255       case COMPARE_CONVERGENCE :  /* Compare values at time t, for 2 iterations */
  256         Val0_Dof = val1[i] ;
  257         Dof_GetRealDofValue(QuantityStorage_S.FunctionSpace->DofData,
  258                 QuantityStorage_S.BasisFunction[0].Dof, &Val1_Dof) ;
  259         break ;
  260       }
  261     }
  262     else
  263       Val0_Dof = Val1_Dof = 0. ;
  264 
  265     val0[i] = Val0_Dof ;  val1[i] = Val1_Dof ;
  266       }  /* for i -> Nbr_Region ... */
  267 
  268 
  269       /*  A n a l y s e   v a l e u r s   . . .  */
  270 
  271       switch (Type_Analyse) {
  272 
  273       case COMPARE_CHANGE :
  274 
  275     switch (ChangeOfState_P->Type) {
  276     case CHANGEOFSTATE_CHANGESIGN :
  277       for (i = 0 ; i < Nbr_Region ; i++) {
  278         if (val0[i] * val1[i] <= 0.)
  279           { *Type_ChangeOfState = CHANGEOFSTATE_CHANGESIGN ;  break ; }
  280       }
  281       break ;
  282 
  283     case CHANGEOFSTATE_CHANGELEVEL :
  284       for (i = 0 ; i < Nbr_Region ; i++) {
  285         if (ChangeOfState_P->Criterion > 0) {
  286           v0 = fabs(val0[i]) ; v1 = fabs(val1[i]) ;
  287           if (((v0 < v1) && (v0*ChangeOfState_P->Criterion < v1)) ||
  288           ((v0 > v1) && (v1*ChangeOfState_P->Criterion < v0)))
  289         { *Type_ChangeOfState = CHANGEOFSTATE_CHANGELEVEL ;  break ; }
  290         }
  291         else {  /* New: Absolute change (Criterion < 0) */
  292           v0 = (val0[i]) ; v1 = (val1[i]) ;
  293           if ( fabs(v1-v0) > fabs(ChangeOfState_P->Criterion) )
  294         { *Type_ChangeOfState = CHANGEOFSTATE_CHANGELEVEL ;  break ; }
  295         }
  296       }  /* Attention: test a affiner ... choix du Criterion ... */
  297       break ;
  298 
  299     case CHANGEOFSTATE_CHANGEREFERENCE :
  300       if (Nbr_Region != 1)
  301         Message::Error("More than 1 Region for ChangeReference not done yet") ;
  302       for (i = 0 ; i < Nbr_Region ; i++) {
  303         if (fabs(val1[i] - val0[i]) >
  304         fabs(ChangeOfState_P->Criterion) *
  305         ((ChangeOfState_P->Criterion > 0.)? fabs(val0[i]) : 1.)) {
  306           *Type_ChangeOfState = ChangeOfState_P->Type ;
  307           *FlagIndex = ChangeOfState_P->FlagIndex ;
  308           if (val1[i] > val0[i])  *FlagIndex *= -1 ;
  309           break ;
  310         }
  311       }
  312       break ;
  313 
  314     case CHANGEOFSTATE_CHANGEREFERENCE2 :
  315       if (Nbr_Region != 1)
  316         Message::Error("More than 1 Region for ChangeReference2 not done yet") ;
  317       for (i = 0 ; i < Nbr_Region ; i++) {
  318         *FlagIndex = ChangeOfState_P->FlagIndex ;
  319         if (val1[i] > val0[i])  *FlagIndex *= -1 ;
  320         if (((struct Expression *)
  321          List_Pointer(Problem_S.Expression, abs(*FlagIndex)))
  322         ->Case.Constant != (*FlagIndex > 0)? 1. : 0. ) {
  323           *Type_ChangeOfState = ChangeOfState_P->Type ;
  324           break ;
  325         }
  326       }
  327       break ;
  328     }
  329 
  330     break ;
  331 
  332       case COMPARE_CONVERGENCE :
  333     Cal_SolutionErrorX(Nbr_Region, val1, val0, &MeanError) ;
  334     if (MeanError > 1.e-8)  *Type_ChangeOfState = !CHANGEOFSTATE_NOCHANGE ;
  335     break ; /* critere a revoir, avant 1.e-14 */
  336 
  337       }
  338 
  339       if (*Type_ChangeOfState != CHANGEOFSTATE_NOCHANGE)  break ;
  340 
  341     }  /* if Nbr_Region > 0 ... */
  342   }  /* for i_COS ... */
  343 
  344   /* Attention: d e b u g  (fprintf)*/
  345 
  346   if ((Type_Analyse == COMPARE_CHANGE &&
  347        (*Type_ChangeOfState != CHANGEOFSTATE_NOCHANGE || !Flag_First)) ||
  348       (Type_Analyse == COMPARE_CONVERGENCE)) {
  349     if (Flag_First) {
  350       for (i = 0 ; i < Nbr_Region ; i++) {
  351     List_Read(Region_L, i, &Num_Region) ; Message::Debug(" %10d",Num_Region) ;
  352       }
  353       for (i = 0 ; i < Nbr_Region ; i++)  Message::Debug(" %.8g", val0[i]) ;
  354     }
  355     for (i = 0 ; i < Nbr_Region ; i++)  Message::Debug(" %.8g", val1[i]) ;
  356     Message::Debug(" t = %.16g, dt = %.16g", Current.Time, Current.DTime) ;
  357     if      (*Type_ChangeOfState == CHANGEOFSTATE_CHANGESIGN)
  358       Message::Debug(" *Sign") ;
  359     else if (*Type_ChangeOfState == CHANGEOFSTATE_CHANGELEVEL)
  360       Message::Debug(" *Level") ;
  361     else if (*Type_ChangeOfState == CHANGEOFSTATE_CHANGEREFERENCE)
  362       Message::Debug(" *Ref (%g %g)",
  363       val0[0]-fabs(ChangeOfState_P->Criterion),
  364       val0[0]+fabs(ChangeOfState_P->Criterion)) ;
  365     else if (*Type_ChangeOfState == CHANGEOFSTATE_CHANGEREFERENCE2)
  366       Message::Debug(" *Ref2 (%g)", val0[0]) ;
  367   }
  368 }
  369 
  370 /* ------------------------------------------------------------------------ */
  371 /*  O p e r a t i o n _ I t e r a t i v e T i m e R e d u c t i o n         */
  372 /* ------------------------------------------------------------------------ */
  373 
  374 void  Operation_IterativeTimeReduction(struct Resolution  * Resolution_P,
  375                        struct Operation   * Operation_P,
  376                        struct DofData     * DofData_P0,
  377                        struct GeoData     * GeoData_P0)
  378 {
  379   int     Num_Iteration, i ;
  380   int     Type_ChangeOfState, Flag_TimeLimLo, Type_LimitHi, FlagIndex ;
  381   double  Time_Previous, DTime0, DTime1 ;
  382   double  Time_LimitLo, DTime_LimitLo, Time_LimitHi ;
  383 
  384   struct Solution   * Solution_P ;
  385   struct Expression * Expression_P ;
  386 
  387 #define TIMELO_OLD   0
  388 #define TIMELO_NEW   1
  389 
  390   Time_Previous = Current.Time - Current.DTime ;
  391   DTime0 = 0. ;  DTime1 = Current.DTime ;
  392 
  393   Flag_TimeLimLo = TIMELO_OLD ;
  394   Time_LimitLo = Time_Previous ;  DTime_LimitLo = Current.DTime ;
  395 
  396   Message::Debug("T I M E   %g (TS #%d, DT %g, Theta %g)",
  397                  Current.Time, (int)Current.TimeStep, Current.DTime, Current.Theta) ;
  398 
  399   Current.SubTimeStep = 0 ;
  400   Treatment_Operation(Resolution_P,
  401               Operation_P->Case.IterativeTimeReduction.Operation,
  402               DofData_P0, GeoData_P0, NULL, NULL) ;
  403   Cal_CompareGlobalQuantity(Operation_P, COMPARE_CHANGE, &Type_ChangeOfState,
  404                 &FlagIndex, 1) ;
  405 
  406   if (Type_ChangeOfState == CHANGEOFSTATE_NOCHANGE) {
  407     Treatment_Operation(Resolution_P,
  408             Operation_P->Case.IterativeTimeReduction.OperationEnd,
  409             DofData_P0, GeoData_P0, NULL, NULL) ;
  410     /* debug */
  411     Cal_CompareGlobalQuantity
  412       (Operation_P, 999, &Type_ChangeOfState, &FlagIndex, 1) ;
  413 
  414     return ;
  415   }
  416 
  417   Time_LimitHi = Current.Time ; /* Sera initialise correctement par apres. */
  418   Type_LimitHi = Type_ChangeOfState ; /* Mais boin, c'est pour la rigueur */
  419 
  420   /* Recherche de l'intervalle de temps [Time_LimitLo, Time_LimitHi] < Criterion
  421      sur lequel un changement d'etat de grandeurs globales specifiees a lieu
  422      (e.g. utilisation pour les circuits avec diodes et thyristors)
  423   */
  424 
  425   for (Num_Iteration = 1 ;
  426        Num_Iteration <= Operation_P->Case.IterativeTimeReduction.NbrMaxIteration ;
  427        Num_Iteration++) {
  428 
  429     if (Type_ChangeOfState == CHANGEOFSTATE_NOCHANGE) {
  430       Flag_TimeLimLo = TIMELO_NEW ;
  431       Time_LimitLo  = Current.Time ;  DTime_LimitLo = Current.DTime ;
  432     }
  433     else {
  434       Time_LimitHi = Current.Time ;
  435       Type_LimitHi = Type_ChangeOfState ;
  436     }
  437 
  438     if (Time_LimitHi - Time_LimitLo <
  439     Operation_P->Case.IterativeTimeReduction.Criterion) {
  440 
  441       if (Type_ChangeOfState != CHANGEOFSTATE_NOCHANGE) {
  442 
  443     if (!(Flag_TimeLimLo == TIMELO_OLD &&
  444           Type_ChangeOfState == CHANGEOFSTATE_CHANGELEVEL) &&
  445         !(Flag_TimeLimLo == TIMELO_OLD)) {
  446       Solution_P = (struct Solution *)
  447         List_Pointer(Current.DofData->Solutions,
  448              List_Nbr(Current.DofData->Solutions)-1) ;
  449       LinAlg_DestroyVector(&Solution_P->x) ;
  450       Free(Solution_P->TimeFunctionValues) ;
  451       Solution_P->SolutionExist = 0 ;
  452       List_Pop(Current.DofData->Solutions) ;  /* Attention: a changer ! */
  453     }
  454 
  455     if (Flag_TimeLimLo == TIMELO_NEW) {
  456       /* Recalcul en Time_LimitLo */
  457       /* Attention: a changer... plutot recuperer solution en Time_LimitLo... */
  458       Message::Debug("==> Re-calculation at Time_LimitLo ... (%.16g)",
  459                          Time_LimitLo) ;
  460       Current.Time  = Time_LimitLo ;  Current.DTime = DTime_LimitLo ;
  461       Current.SubTimeStep++ ;
  462 
  463       Treatment_Operation(Resolution_P,
  464                   Operation_P->Case.IterativeTimeReduction.Operation,
  465                   DofData_P0, GeoData_P0, NULL, NULL) ;
  466     }
  467       }
  468 
  469       if (Flag_TimeLimLo == TIMELO_NEW ||
  470       (Flag_TimeLimLo == TIMELO_OLD &&
  471        Type_ChangeOfState == CHANGEOFSTATE_CHANGELEVEL)) {
  472     Treatment_Operation(Resolution_P,
  473                 Operation_P->Case.IterativeTimeReduction.OperationEnd,
  474                 DofData_P0, GeoData_P0, NULL, NULL) ;
  475     /* debug */
  476     Cal_CompareGlobalQuantity
  477       (Operation_P, 999, &Type_ChangeOfState, &FlagIndex, 1) ;
  478       }
  479 
  480       if (Type_LimitHi == CHANGEOFSTATE_CHANGESIGN ||
  481       Type_LimitHi == CHANGEOFSTATE_CHANGEREFERENCE ||
  482       Type_LimitHi == CHANGEOFSTATE_CHANGEREFERENCE2) {
  483 
  484     if (Flag_TimeLimLo == TIMELO_NEW)  Current.TimeStep += 1. ;
  485 
  486     if (Type_LimitHi == CHANGEOFSTATE_CHANGEREFERENCE ||
  487         Type_LimitHi == CHANGEOFSTATE_CHANGEREFERENCE2) {
  488       Expression_P = (struct Expression *)
  489         List_Pointer(Problem_S.Expression, abs(FlagIndex)) ;
  490       Expression_P->Case.Constant = (FlagIndex > 0)? 1. : 0. ;
  491       /*
  492         Expression_P->Case.Constant =
  493         (double)(!((int)Expression_P->Case.Constant)) ;
  494       */
  495       Message::Debug("===> Flag -> %g", Expression_P->Case.Constant) ;
  496     }
  497 
  498     if (Operation_P->Case.IterativeTimeReduction.Flag)
  499       Current.Theta = 1. ;
  500     /* New: Theta is also changed for this time !
  501        OK because dt is then very small also in this case ! */
  502     Current.Time  = Time_LimitHi ;
  503     Current.DTime = Time_LimitHi - Time_LimitLo ;
  504     Current.SubTimeStep++ ;
  505 
  506     Message::Debug("==> iterations for TimeHi ...") ;
  507 
  508     i = 0 ;
  509     do {
  510       i ++ ;
  511       Treatment_Operation(Resolution_P,
  512                   Operation_P->Case.IterativeTimeReduction.Operation,
  513                   DofData_P0, GeoData_P0, NULL, NULL) ;
  514       Cal_CompareGlobalQuantity
  515         (Operation_P, COMPARE_CONVERGENCE, &Type_ChangeOfState, &FlagIndex, 0) ;
  516     } while ((Flag_TimeLimLo == TIMELO_NEW && i == 1) ||
  517          (Type_ChangeOfState != CHANGEOFSTATE_NOCHANGE && i < 9)) ;
  518     /* Attention: critere (NbrMax 9) a revoir */
  519 
  520     Treatment_Operation(Resolution_P,
  521                 Operation_P->Case.IterativeTimeReduction.OperationEnd,
  522                 DofData_P0, GeoData_P0, NULL, NULL) ;
  523     /* debug */
  524     Cal_CompareGlobalQuantity
  525       (Operation_P, 999, &Type_ChangeOfState, &FlagIndex, 1) ;
  526 
  527     if (Operation_P->Case.IterativeTimeReduction.Flag) { /* Attention: Test */
  528       Message::Debug("=====> Theta = %g -> 1.", Current.Theta) ;
  529       Flag_NextThetaFixed = 1 ;
  530       Current.Theta = 1. ;
  531       if (Operation_P->Case.IterativeTimeReduction.Flag > 0) {
  532         Current.DTime *= (double)Operation_P->Case.IterativeTimeReduction.Flag ;
  533         Flag_NextThetaFixed = 2 ;  /* Theta is fixed, DTime is also fixed */
  534       }
  535     }
  536       }
  537 
  538       break ;  /* Out of loop 'for Num_Iteration' */
  539     }  /* if Time_LimitHi - Time_LimitLo << ... */
  540 
  541 
  542     if (Operation_P->Case.IterativeTimeReduction.DivisionCoefficient > 0.) {
  543       if (Type_ChangeOfState == CHANGEOFSTATE_NOCHANGE)  DTime0 += DTime1 ;
  544       DTime1 /= Operation_P->Case.IterativeTimeReduction.DivisionCoefficient ;
  545     }
  546     else { /* Technique de Pkp ... "un peu trop prudente" */
  547       if (Type_ChangeOfState == CHANGEOFSTATE_NOCHANGE)  DTime0 += DTime1 ;
  548       else  DTime1 /=
  549           fabs(Operation_P->Case.IterativeTimeReduction.DivisionCoefficient) ;
  550     }
  551 
  552     Current.DTime = DTime0 + DTime1 ;
  553     Current.Time = Time_Previous + Current.DTime ;
  554     Current.SubTimeStep++ ;
  555 
  556     Solution_P = (struct Solution *)
  557       List_Pointer(Current.DofData->Solutions,
  558            List_Nbr(Current.DofData->Solutions)-1) ;
  559     LinAlg_DestroyVector(&Solution_P->x) ;
  560     Free(Solution_P->TimeFunctionValues) ;
  561     Solution_P->SolutionExist = 0 ;
  562     List_Pop(Current.DofData->Solutions) ;  /* Attention: a changer ! */
  563 
  564     Treatment_Operation(Resolution_P,
  565             Operation_P->Case.IterativeTimeReduction.Operation,
  566             DofData_P0, GeoData_P0, NULL, NULL) ;
  567     Cal_CompareGlobalQuantity(Operation_P, COMPARE_CHANGE, &Type_ChangeOfState,
  568                   &FlagIndex, 0) ;
  569   }  /* for Num_Iteration ... */
  570 }