~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/ode/Homotopy.cpp

• browse files at revision 1693
• compare with another revision
• download diff
• download tarball
• view history from revision 1693

Show diffs side-by-side

added

removed

src/kernel/ode/Homotopy.cpp

1

 

// Copyright (C) 2005 Anders Logg.

 

1

// Copyright (C) 2005-2006 Anders Logg.

2

2

// Licensed under the GNU GPL Version 2.

3

3

//

 

4

// Modified by Garth N. Wells, 2006.

 

5

//

4

6

// First added:  2005

5

 

// Last changed: 2005-12-20

 

7

// Last changed: 2006-08-22

6

8

 

7

9

#include <stdio.h>

 

10

#include <limits>

8

11

#include <dolfin/dolfin_log.h>

9

12

#include <dolfin/dolfin_math.h>

10

13

#include <dolfin/ParameterSystem.h>

11

 

#include <dolfin/LU.h>

12

 

#include <dolfin/GMRES.h>

13

14

#include <dolfin/ComplexODE.h>

14

15

#include <dolfin/HomotopyJacobian.h>

15

16

#include <dolfin/HomotopyODE.h>

20

21

//-----------------------------------------------------------------------------

21

22

Homotopy::Homotopy(uint n)

22

23

  : tol(0), n(n), M(0), maxiter(0), maxpaths(0), maxdegree(0),

23

 

    divtol(0), monitor(false), random(false), solver(0), filename(""),

24

 

    mi(0), ci(0), tmp(0), x(2*n),

 

24

    divtol(0), monitor(false), random(false), 

 

25

    filename(""), mi(0), ci(0), tmp(0), x(2*n),

25

26

    degree_adjusted("Adjusting degree of equation, maximum reached.")

26

27

{

27

28

  dolfin_info("Creating homotopy for system of size %d.", n);

28

29

  

29

30

  // We should not solve the dual problem

30

 

  set("solve dual problem", false);

 

31

  set("ODE solve dual problem", false);

31

32

 

32

33

  // System is implicit

33

 

  set("implicit", true);

 

34

  set("ODE implicit", true);

34

35

 

35

36

  // Get divergence tolerance

36

37

  divtol = get("homotopy divergence tolerance");

44

45

    dolfin_info("Using random initial system for homotopy.");

45

46

 

46

47

  // Get maximum number of iterations

47

 

  maxiter = get("maximum iterations");

 

48

  maxiter = get("ODE maximum iterations");

48

49

 

49

50

  // Get maximum number of paths

50

51

  maxpaths = get("homotopy maximum size");

62

63

  // FIXME: Maybe this should be a parameter?

63

64

  tol = get("homotopy solution tolerance");

64

65

  

65

 

  // Choose solver

66

 

  //solver = new GMRES();

67

 

  //((GMRES *) solver)->setAtol(0.1*tol);

68

 

  //((GMRES *) solver)->setReport(false);

69

 

  solver = new LU();

70

 

 

71

66

  // Initialize array mi

72

67

  mi = new uint[n];

73

68

  for (uint i = 0; i < n; i++)

87

82

  if ( mi ) delete [] mi;

88

83

  if ( ci ) delete [] ci;

89

84

  if ( tmp ) delete [] tmp;

90

 

  if ( solver ) delete solver;

91

85

  

92

86

  for (unsigned int i = 0; i < zs.size(); i++)

93

87

    delete [] zs[i];

109

103

    dolfin_info("\nComputing path number %d out of %d.", m + 1, M);

110

104

 

111

105

    // Change name of output file for each path

112

 

    sprintf(filename, "primal_%d.m", m);

113

 

    set("ode solution file name", filename);

 

106

    sprintf(filename, "solution_%u.py", m);

 

107

    set("ODE solution file name", filename);

114

108

 

115

109

    // Compute the component paths from global path number

116

110

    computePath(m);

151

145

  return zs;

152

146

}

153

147

//-----------------------------------------------------------------------------

154

 

complex Homotopy::z0(uint i)

 

148

void Homotopy::z0(complex z[])

155

149

{

156

 

  const real pp = static_cast<real>(adjustedDegree(i));

157

 

  const real mm = static_cast<real>(mi[i]);

158

 

  const complex c = ci[i];

159

 

  

160

 

  // Pick root number m of equation z_i^(p + 1) = c_i

161

 

  real r = std::pow(std::abs(c), 1.0/(pp + 1.0));

162

 

  real a = std::arg(c) / (pp + 1.0);

163

 

  complex z = std::polar(r, a + mm/(pp + 1.0)*2.0*DOLFIN_PI);

164

 

  

165

 

  return z;

 

150

  for (uint i = 0; i < n; i++)

 

151

  {

 

152

    const real pp = static_cast<real>(adjustedDegree(i));

 

153

    const real mm = static_cast<real>(mi[i]);

 

154

    const complex c = ci[i];

 

155

    

 

156

    // Pick root number m of equation z_i^(p + 1) = c_i

 

157

    real r = std::pow(std::abs(c), 1.0/(pp + 1.0));

 

158

    real a = std::arg(c) / (pp + 1.0);

 

159

    z[i] = std::polar(r, a + mm/(pp + 1.0)*2.0*DOLFIN_PI);

 

160

  }

166

161

}

167

162

//-----------------------------------------------------------------------------

168

163

void Homotopy::G(const complex z[], complex y[])

273

268

bool Homotopy::computeSolution(HomotopyODE& ode)

274

269

{

275

270

  // Create right-hand side and increment vector

276

 

  Vector F(2*n), dx(2*n);

 

271

  uBlasVector F(2*n), dx(2*n);

277

272

 

278

273

  // Create matrix-free Jacobian

279

274

  HomotopyJacobian J(ode, x);

280

 

  J.init(dx, dx);

281

275

  

282

276

  cout << "Starting point:     x = ";

283

277

  x.disp();

289

283

    feval(F, ode);

290

284

 

291

285

    // Check convergence

292

 

    real r = F.norm(Vector::linf);

 

286

    real r = F.norm(uBlasVector::linf);

293

287

    //cout << "r = " << r << ": x = "; x.disp();

294

288

    if ( r < tol )

295

289

    {

297

291

      x.disp();

298

292

      return true;

299

293

    }

300

 

    

301

 

    //cout << "x = "; x.disp();

302

 

    //cout << "F = "; F.disp();

303

 

    //cout << "dx = "; dx.disp();

304

 

    //cout << "J = "; J.disp(false);

305

 

    //cout << endl;

306

294

 

307

 

    // Solve linear system, seems like we need to scale the right-hand

308

 

    // side to make it work with the PETSc GMRES solver

309

 

    //r += DOLFIN_EPS;

310

 

    //F /= r;

311

 

    //solver->solve(J, dx, F);

312

 

    //dx *= r;

313

 

    

314

 

    // Solve linear system using LU factorization

315

 

    solver->solve(J, dx, F);

 

295

    // Solve linear system

 

296

    solver.solve(J, dx, F);

316

297

 

317

298

    // Subtract increment

318

299

    x -= dx;

319

 

    //x.disp();

320

300

  }

321

301

 

322

302

  dolfin_warning("Solution did not converge.");

326

306

void Homotopy::saveSolution()

327

307

{

328

308

  // Copy values to complex array

329

 

  real* xx = x.array();

330

309

  complex* z = new complex[n];

331

310

  for (uint i = 0; i < n; i++)

332

 

    z[i] = complex(xx[2*i], xx[2*i + 1]);

333

 

  x.restore(xx);

 

311

    z[i] = complex(x(2*i), x(2*i + 1));

334

312

 

335

313

  // Allow user to modify solution

336

314

  modify(z);

387

365

  fclose(fp);

388

366

}

389

367

//-----------------------------------------------------------------------------

390

 

void Homotopy::feval(Vector& F, ComplexODE& ode)

 

368

void Homotopy::feval(uBlasVector& F, ComplexODE& ode)

391

369

{

392

370

  // Reuse the right-hand side of the ODE so we don't have to reimplement

393

371

  // the mapping from complex to real numbers

394

372

 

395

 

  // Get arrays

396

 

  const real* xx = x.array();

397

 

  real* FF = F.array();

398

 

  

399

373

  // Evaluate F at current x

400

 

  ode.f(xx, 0.0, FF);

401

 

 

402

 

  // Restore arrays

403

 

  x.restore(xx);

404

 

  F.restore(FF);

 

374

  ode.f(x, 0.0, F);

405

375

}

406

376

//-----------------------------------------------------------------------------

• Older »

Loggerhead is a web-based interface for Breezy
Version: 2.0.1