OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
linesearch.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "linesearch.h"
36 #include "timestep.h"
37 #include "floatarray.h"
38 #include "intarray.h"
39 #include "mathfem.h"
40 #include "nmstatus.h"
41 #include "engngm.h"
42 
43 namespace oofem {
45  NumericalMethod(d, m)
46 {
47  max_iter = 10;
48  ls_tolerance = 0.80;
49  amplifFactor = 2.5;
50  maxEta = 4.0;
51  minEta = 0.2;
52 }
53 
56  IntArray &eqnmask, double lambda, double &etaValue, LS_status &status, TimeStep *tStep)
57 {
58  int ico, ils, neq = r.giveSize();
59  double s0;
60 
61  FloatArray g(neq), rb(neq);
62  // Compute inner product at start and stop if positive
63  g = R;
64  g.times(lambda);
65  if ( R0 ) {
66  g.add(* R0);
67  }
68 
69  g.subtract(F);
70 
71  for ( auto eq : eqnmask ) {
72  g.at( eq ) = 0.0;
73  }
74 
75  s0 = ( -1.0 ) * g.dotProduct(dr);
76  if ( s0 >= 0.0 ) {
77  //printf ("\nLineSearchNM::solve starting inner product uphill, val=%e",s0);
78  OOFEM_LOG_DEBUG("LS: product uphill, eta=%e\n", 1.0);
79  r.add(dr);
80  tStep->incrementStateCounter(); // update solution state counter
82  etaValue = 1.0;
83  status = ls_ok;
84  return NM_Success;
85  }
86 
87  // keep original total displacement r
88  rb = r;
89 
90  eta.resize(this->max_iter + 1);
91  prod.resize(this->max_iter + 1);
92  // prepare starting product ratios and step lengths
93  prod.at(1) = 1.0;
94  eta.at(1) = 0.0;
95  eta.at(2) = 1.0;
96  // following counter shows how many times the max or min step length has been reached
97  ico = 0;
98 
99  // begin line search loop
100  for ( ils = 2; ils <= this->max_iter; ils++ ) {
101  // update displacements
102  r = rb;
103  r.add(this->eta.at(ils), dr);
104 
105  tStep->incrementStateCounter(); // update solution state counter
106  // update internal forces according to new state
108  // compute out-of balance forces g in new state
109  g = R;
110  g.times(lambda);
111  if ( R0 ) {
112  g.add(* R0);
113  }
114 
115  g.subtract(F);
116 
117  for ( auto eq : eqnmask ) {
118  g.at( eq ) = 0.0;
119  }
120 
121  // compute current inner-product ratio
122  double si = ( -1.0 ) * g.dotProduct(dr) / s0;
123  prod.at(ils) = si;
124 
125  // check if line-search tolerance is satisfied
126  if ( fabs(si) < ls_tolerance ) {
127  dr.times( this->eta.at(ils) );
128  //printf ("\nLineSearchNM::solve tolerance satisfied for eta=%e, ils=%d", eta.at(ils),ils);
129  OOFEM_LOG_DEBUG( "LS: ils=%d, eta=%e\n", ils, eta.at(ils) );
130 
131  etaValue = eta.at(ils);
132  status = ls_ok;
133  return NM_Success;
134  }
135 
136  // call line-search routine to get new estimate of eta.at(ils)
137  this->search(ils, prod, eta, this->amplifFactor, this->maxEta, this->minEta, ico);
138  if ( ico == 2 ) {
139  break; // exit the loop
140  }
141  } // end line search loop
142 
143  // exceeded no of ls iterations of ls failed
144  //if (ico == 2) printf("\nLineSearchNM::solve max or min step length has been reached two times");
145  //else printf("\nLineSearchNM::solve reached max number of ls searches");
146  OOFEM_LOG_DEBUG( "LS: ils=%d, ico=%d, eta=%e\n", ils, ico, eta.at(ils) );
147  /* update F before */
148  r = rb;
149  r.add(dr);
150 
151  tStep->incrementStateCounter(); // update solution state counter
153  etaValue = 1.0;
154  status = ls_failed;
155  return NM_NoSuccess;
156 }
157 
158 
159 void
161  double maxetalim, double minetalim, int &ico)
162 {
163  int ineg = 0;
164  double etaneg = 1.0;
165  double etamax = 0.0;
166 
167 
168  // obtain ineg (number of previous line search iteration with negative ratio nearest to origin)
169  // as well as max previous step length, etamax
170 
171  for ( int i = 1; i <= istep; i++ ) {
172  etamax = max( etamax, eta.at(i) );
173  if ( prod.at(i) >= 0.0 ) {
174  continue;
175  }
176 
177  if ( eta.at(i) >= etaneg ) {
178  continue;
179  }
180 
181  etaneg = eta.at(i);
182  ineg = i;
183  }
184 
185  if ( ineg ) {
186  // allow interpolation
187  // first find ipos (position of previous s-l with positive ratio that is
188  // closest to ineg (but with smaller s-l)
189  int ipos = 1;
190  for ( int i = 1; i <= istep; i++ ) {
191  if ( prod.at(i) <= 0.0 ) {
192  continue;
193  }
194 
195  if ( eta.at(i) > eta.at(ineg) ) {
196  continue;
197  }
198 
199  if ( eta.at(i) < eta.at(ipos) ) {
200  continue;
201  }
202 
203  ipos = i;
204  }
205 
206  // interpolate to get step-length
207  double etaint = ( prod.at(ineg) * eta.at(ipos) - prod.at(ipos) * eta.at(ineg) ) / ( prod.at(ineg) - prod.at(ipos) );
208  // alternativelly get eta ensuring reasonable change
209  double etaalt = eta.at(ipos) + 0.2 * ( eta.at(ineg) - eta.at(ipos) );
210  etaint = max(etaint, etaalt);
211  if ( etaint < minetalim ) {
212  etaint = minetalim;
213  if ( ico == 1 ) {
214  ico = 2;
215  } else {
216  ico = 1;
217  }
218  }
219 
220  eta.at(istep + 1) = etaint;
221  return;
222  } else { // ineq == 0
223  // allow extrapolation
224  double etamaxstep = amp * etamax;
225  // extrapolate between current and previous
226  double etaextrap = ( prod.at(istep) * eta.at(istep - 1) - prod.at(istep - 1) * eta.at(istep) ) /
227  ( prod.at(istep) - prod.at(istep - 1) );
228  eta.at(istep + 1) = etaextrap;
229  // check if in limits
230  if ( ( etaextrap <= 0.0 ) || ( etaextrap > etamaxstep ) ) {
231  eta.at(istep + 1) = etamaxstep;
232  }
233 
234  if ( ( eta.at(istep + 1) > maxetalim ) && ( ico == 1 ) ) {
235  ico = 2;
236  return;
237  }
238 
239  if ( ( eta.at(istep + 1) > maxetalim ) ) {
240  ico = 1;
241  eta.at(istep + 1) = maxetalim;
242  }
243  }
244 }
245 
248 {
249  IRResultType result; // Required by IR_GIVE_FIELD macro
250 
251  /* default values set in constructor
252  * ls_tolerance = 0.80;
253  * amplifFactor = 2.5;
254  * maxEta = 4.0;
255  */
257  if ( ls_tolerance < 0.6 ) {
258  ls_tolerance = 0.6;
259  }
260 
261  if ( ls_tolerance > 0.95 ) {
262  ls_tolerance = 0.95;
263  }
264 
266  if ( amplifFactor < 1.0 ) {
267  amplifFactor = 1.0;
268  }
269 
270  if ( amplifFactor > 10.0 ) {
271  amplifFactor = 10.0;
272  }
273 
275  if ( maxEta < 1.5 ) {
276  maxEta = 1.5;
277  }
278 
279  if ( maxEta > 15.0 ) {
280  maxEta = 15.0;
281  }
282 
283  //printf ("\nLineSearchNM::initializeFrom: tol=%e, ampl=%e, maxEta=%e\n",
284  // ls_tolerance, amplifFactor,maxEta);
285 
286  return IRRT_OK;
287 }
288 } // end namespace oofem
FloatArray prod
Definition: linesearch.h:70
#define NM_Success
Numerical method exited with success.
Definition: nmstatus.h:47
Class and object Domain.
Definition: domain.h:115
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
This base class is an abstraction for numerical algorithm.
Definition: nummet.h:80
void incrementStateCounter()
Updates solution state counter.
Definition: timestep.h:190
unsigned long NM_Status
Mask defining NumMetod Status; which can be asked after finishing computation by Numerical Method...
Definition: nmstatus.h:44
#define OOFEM_LOG_DEBUG(...)
Definition: logger.h:128
Class implementing an array of integers.
Definition: intarray.h:61
LineSearchNM(Domain *d, EngngModel *m)
Constructor.
Definition: linesearch.C:44
#define NM_NoSuccess
Numerical method failed to solve problem.
Definition: nmstatus.h:48
void search(int istep, FloatArray &prod, FloatArray &eta, double amp, double maxeta, double mineta, int &status)
Definition: linesearch.C:160
Domain * domain
Pointer to domain.
Definition: nummet.h:84
virtual void updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
Updates components mapped to numerical method if necessary during solution process.
Definition: engngm.C:1485
Class representing vector of real numbers.
Definition: floatarray.h:82
#define _IFT_LineSearchNM_lsearchamp
Definition: linesearch.h:46
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
virtual NM_Status solve(FloatArray &r, FloatArray &dr, FloatArray &F, FloatArray &R, FloatArray *R0, IntArray &eqnmask, double lambda, double &etaValue, LS_status &status, TimeStep *tStep)
Solves the line search optimization problem in the form of , The aim is to find so that the has dec...
Definition: linesearch.C:55
virtual IRResultType initializeFrom(InputRecord *ir)
Definition: linesearch.C:247
Class representing the general Input Record.
Definition: inputrecord.h:101
FloatArray eta
Definition: linesearch.h:69
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
#define _IFT_LineSearchNM_lsearchtol
Definition: linesearch.h:45
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
EngngModel * engngModel
Pointer to engineering model.
Definition: nummet.h:86
#define _IFT_LineSearchNM_lsearchmaxeta
Definition: linesearch.h:47
Class representing solution step.
Definition: timestep.h:80
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011