Mercurialsp / file revision
summary | shortlog | changelog | graph | tags | files | changeset | file | revisions | annotate | diff | raw
src/test_chafe.cpp
author Alexey Ozeritsky <aozeritsky@gmail.com>
Fri Jun 11 21:50:18 2010 +0400 (20 months ago)
changeset 142 2aa3b89e9ccc
parent 111993c76371d3e
child 176acf2a0db4441
permissions -rw-r--r--
fix 
     1 #include <math.h>

     2 #include <stdlib.h>

     3 #include <stdio.h>

     4 

     5 #include <vector>

     6 

     7 #include "chafe.h"

     8 #include "linal.h"

     9 

    10 #ifdef max

    11 #undef max

    12 #endif

    13 

    14 using namespace std;

    15 using namespace linal;

    16 

    17 static inline double max (double a, double b)

    18 {

    19 	return (a > b) ? a : b;

    20 }

    21 

    22 double ans (double x, double y, double t)

    23 {

    24 	return x*sin (y + t) *ipow (cos (x), 4);

    25 }

    26 

    27 double f (double x, double y, double t,

    28           double mu, double sigma)

    29 {

    30 	return ipow (cos (x), 2) *

    31 	       (x*cos (y + t) *ipow (cos (x), 2)

    32 	        + sigma*sin (y + t) *ipow (cos (x), 2) *x + 9*mu*sin (y + t) *sin (x) *cos (x)

    33 	        - 15*mu*sin (y + t) *x + 20*mu*sin (y + t) *x*ipow (cos (x), 2) );		

    34 }

    35 

    36 void solve()

    37 {

    38 	long nlat = 19;

    39 	long nlon = 36;

    40 

    41 	SphereChafeConf conf;

    42 	conf.nlat  = nlat;

    43 	conf.nlon  = nlon;

    44 	conf.mu    = 1.0;

    45 	conf.sigma = +70.0;

    46 	conf.tau   = 0.01;

    47 	conf.theta = 0.5;

    48 	conf.rp    = f;

    49 

    50 	double dlat = M_PI / (nlat - 1);

    51 	double dlon = 2. * M_PI / nlon;

    52 	double t = 0;

    53 

    54 	int i, j, it = 0;

    55 

    56 	vector < double > u (nlat * nlon);

    57 	vector < double > v (nlat * nlon);

    58 	vector < double > r (nlat * nlon);

    59 

    60 	SphereChafe chafe (conf);

    61 

    62 	double nev1 = 0;

    63 

    64 	for (i = 0; i < nlat; ++i)

    65 	{

    66 		for (j = 0; j < nlon; ++j)

    67 		{

    68 			double phi    = -0.5 * M_PI + i * dlat;

    69 			double lambda = j * dlon;

    70 

    71 			r[i * nlon + j] = ans (phi, lambda, t);

    72 		}

    73 	}

    74 

    75 	while (true)

    76 	{

    77 		chafe.calc (&u[0], &r[0], t);

    78 		t += conf.tau;

    79 

    80 		if (it % 100 == 0)

    81 		{

    82 			nev1 = 0.0;

    83 			for (i = 0; i < nlat; ++i)

    84 			{

    85 				for (j = 0; j < nlon; ++j)

    86 				{

    87 					double phi    = -0.5 * M_PI + i * dlat;

    88 					double lambda = j * dlon;

    89 

    90 					v[i * nlon + j] = ans (phi, lambda, t);

    91 					nev1 = max (nev1, fabs (u[i * nlon + j] - v[i * nlon + j] ) );

    92 

    93 //					fprintf(stderr, "%.16lf %.16lf \n", u[i * nlon + j], v[i * nlon + j]);

    94 				}

    95 			}

    96 

    97 			fprintf (stderr, "nev1=%.16le \n", nev1);

    98 		}

    99 		r.swap(u);

   100 

   101 		it += 1;

   102 	}

   103 }

   104 

   105 int main (int argc, char * argv[])

   106 {

   107 	solve();

   108 }

   109
sp