← Back to branch summary

~andy-terrel/junk/dolfwave

« back to all changes in this revision

Viewing changes to tests/1HD1S/solitaryKdVcdgP1P2/main.cpp

  • Committer: Nuno Lopes
  • Date: 2013-01-04 20:33:34 UTC
  • Revision ID: ndlopes@gmail.com-20130104203334-fi37ndqxbhrn8jib
Snapshot. Updated Changelog

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/1HD1S/solitaryKdVcdgP1P2/main.cpp
1
 
// Copyright (C) 2012 Nuno David Lopes.-----------------------------------
2
 
// e-mail: ndl@ptmat.fc.ul.pt
3
 
// http://ptmat.fc.ul.pt/~ndl/
4
 
// Licensed under the GNU LGPL Version 3.0 or later.
5
 
// This program provides a demo of the DOLFWAVE library.
6
 
// Solitary wave evolution by adimensional generic
7
 
// KdV-BBM with CDG method
8
 
 
9
 
#include <dolfwave.h>
10
 
using namespace dolfin;
11
 
using namespace dolfin::dolfwave;
12
 
 
13
 
/* Global Variables*/
14
 
const double A=-20.; //
15
 
const double B=200.;  // $x\in[A.B]$
16
 
const double eps=0.1;
17
 
const double mu=0.31623;
18
 
 
19
 
class ElevationInitKdV: public Expression
20
 
{
21
 
public:
22
 
  ElevationInitKdV(const double& delta0,const double& epsilon,
23
 
                   const double& delta1,const double& delta2):
24
 
    delta0(delta0), epsilon(epsilon), delta1(delta1), delta2(delta2) {}
25
 
  void eval(Array<double> & values,const Array<double> & x) const
26
 
  {
27
 
    /*
28
 
       Definition of the solitary wave solution
29
 
       for the adimensional generic KdV-BBM equations
30
 
    */
31
 
 
32
 
    double c=delta0+delta0/80.; //Free parameter for speed/wave height
33
 
    double b=3./2.*epsilon; //Nonlinear coeff.
34
 
    double Psi=3.*(c-delta0)/b;
35
 
    double Theta=1./2.*sqrt((c-delta0)/(delta2-c*delta1));
36
 
    values[0]=Psi/sqr(cosh(Theta*x[0]));
37
 
  }
38
 
private:
39
 
  const double &delta0;
40
 
  const double &epsilon;
41
 
  const double &delta1;
42
 
  const double &delta2;
43
 
};
44
 
// --------------------------------------------------------------------
45
 
 
46
 
int main(int argc, char* argv[])
47
 
{
48
 
  bool swtch=true;
49
 
  if (argc!=4)
50
 
    {
51
 
      info("Wrong number of arguments");
52
 
      info("argv[1]='2' for KdV-BBM or argv[1]='1' for KdV");
53
 
      info("argv[2]=penalty_tau");
54
 
      info("argv[3]=output_dir");
55
 
      exit(EXIT_FAILURE);
56
 
    }
57
 
  std::stringstream output_dir;
58
 
  output_dir<<argv[3];
59
 
  //  Options, Initialization, Mesh and User functions
60
 
  Dolfwave dw(120000 /*MaxSteps*/,
61
 
              0.001 /*TimeStep*/,
62
 
              1000 /*WriteGap*/,
63
 
              "aKdVcdgP1P2_1D" /*FormName*/,
64
 
              "LU_P" /* uBLAS AlgBackEND*/,
65
 
              "gnuplot" /*Plotter*/,
66
 
              output_dir.str() /*OutDir*/,
67
 
              "pvd"/*OutFormat*/);
68
 
 
69
 
  if (atoi(argv[1])==1)
70
 
    dw.SetParameter("KdVPade",kdv_21_pade);
71
 
  else if (atoi(argv[1])==2)
72
 
    dw.SetParameter("KdVPade",kdv_22_pade);
73
 
 
74
 
  double kdvp=parameters("dolfwave_parameters")["KdVPade"];
75
 
  double delta1=(kdvp-1.)*mu*mu/6.; //Sec.-ord. disp. coeff.
76
 
  double delta2=kdvp*mu*mu/6.; //Third-ord. disp. coeff.
77
 
 
78
 
  dw.SetParameter("aKdVDelta0",1.0);
79
 
  dw.SetParameter("aKdVEpsilon",eps);
80
 
  dw.SetParameter("aKdVDelta1",delta1);
81
 
  dw.SetParameter("aKdVDelta2",delta2);
82
 
  dw.SetParameter("PenaltyTau",atof(argv[2]));
83
 
  info(parameters,true);// Some info on the used parameters
84
 
  Interval mesh(1001,A,B);
85
 
  ElevationInitKdV eta_init(parameters("dolfwave_parameters")["aKdVDelta0"],
86
 
                            parameters("dolfwave_parameters")["aKdVEpsilon"],
87
 
                            parameters("dolfwave_parameters")["aKdVDelta1"],
88
 
                            parameters("dolfwave_parameters")["aKdVDelta2"]);
89
 
  //initial condition for elevation
90
 
  dw.SpacesFunctionsVectorsInit(mesh); // Init. Spaces Functions & Vectors
91
 
  dw.BilinearFormInit(mesh);//BilinearForm a initialization
92
 
  dw.MatricesAssemble( ); //Matrices M and Mf initialization
93
 
  dw.LinearFormsInit( );
94
 
  dw.InitialCondition(eta_init); //Initial Conditions L2 projection
95
 
  Constant zero(0.0);
96
 
  dw.Volume(mesh,zero);
97
 
  dw.LUFactorization(true ); //Re-use factorization
98
 
  dw.SaveEta(0);
99
 
  dw.PreviewEta(0);
100
 
  dw.RKInit("exp4"); // Explicit 4th-order Runge-Kutta
101
 
  dw.RKSolve();//Runge-Kutta for 3 initial steps
102
 
  dw.PCInit(mesh,true);
103
 
 
104
 
  for(int i=4; i<dw.MaxSteps+1;i++)
105
 
    {
106
 
      dw.PCSolve( );
107
 
      if (!(i%dw.WriteGap))
108
 
        { // Plots
109
 
          dw.PreviewEta();
110
 
          dw.SaveEta();
111
 
          dw.Volume(mesh,zero);
112
 
          dw.Info(i,true);// Output some extra information on memory usage
113
 
        }
114
 
      dw.Info(i,false);
115
 
    }
116
 
  return (EXIT_SUCCESS);
117
 
}