GALA2D/Fields.cpp at master · pseudospectators/GALA2D · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
/*************************************************************************************************************************
 *
 * Class Solver
 *
 * Functions for initial condition and velocity
 *
 *************************************************************************************************************************/

#define _USE_MATH_DEFINES

#include <iostream>
#include <fstream>
#include <cassert>
#include <cmath>
#include <ctime>
#include "QuadTree.h"
#include "QuadNode.h"
#include "Point.h"
#include "Solver.h"
using namespace std;

/*-------------------------------------------------------------------------*\
\*-------------------------------------------------------------------------*/

// Initial condition
// For given x and y (correspond to x_1 and x_2 in [1])
// computes initial condition phi (corresponds to u in [1]),
// its first derivatives and its cross-derivative
void Solver::inco(double x, double y, double& phi, double& phix, double& phiy, double& phixy)
{

   // Periodized Gaussian, eq (3.2) in [1]
   // Initialize to zero
   double phi0 = 0.0;
   double phi0x = 0.0;
   double phi0y = 0.0;
   double phi0xy = 0.0;
   // Sum up Gaussians
   for (int ix=-30; ix<=30; ix++) {
      for (int iy=-30; iy<=30; iy++) {
         double xi = x - double(ix);
         double yi = y - double(iy);
         phi0 = phi0 + exp( -10.0*((xi-0.5)*(xi-0.5)+(yi-0.75)*(yi-0.75)) );
         phi0x = phi0x - 20.0*(xi-0.5) * exp( -10.0*((xi-0.5)*(xi-0.5)+(yi-0.75)*(yi-0.75)) );
         phi0y = phi0y - 20.0*(yi-0.75) * exp( -10.0*((xi-0.5)*(xi-0.5)+(yi-0.75)*(yi-0.75)) );
         phi0xy = phi0xy + 400.0*(yi-0.75)*(xi-0.5) * exp( -10.0*((xi-0.5)*(xi-0.5)+(yi-0.75)*(yi-0.75)) );
   } }
   // Return values
   phi = phi0;
   phix = phi0x;
   phiy = phi0y;
   phixy = phi0xy;

}

/*-------------------------------------------------------------------------*\
\*-------------------------------------------------------------------------*/

// Advection velocity
// For given x, y and time t
// compute the velocity components vx and vy (correspond to a_1 and a_2 in [1])
void Solver::velocity(double x, double y, double t, double& vx, double& vy)
{
   // Swirl test, eq (3.3) in [1]
   double T = 10.0;
   double xi = x;
   double yi = y;
   // TODO: M_PI may need to be redefined for double precision
   vx = cos(M_PI*t/T) * sin(M_PI*xi) * sin(M_PI*xi) * sin(2*M_PI*yi);
   vy = - cos(M_PI*t/T) * sin(2*M_PI*xi) * sin(M_PI*yi) * sin(M_PI*yi);
}

/*-------------------------------------------------------------------------*\
\*-------------------------------------------------------------------------*/