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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include "ectrans/transi.h"
#include "transi_test.h"
// ----------------------------------------------------------------------------
double randomDouble() {
uint64_t r53 = ((uint64_t)(rand()) << 21) ^ (rand() >> 2);
return (double)r53 / 9007199254740991.0; // 2^53 - 1
}
// ----------------------------------------------------------------------------
void test_lam_dirtrans_adjoint() {
double adjoint_tol = 1.e-6;
printf("test_lam_dirtrans_adjoint()\n");
struct Trans_t trans;
TRANS_CHECK( trans_new(&trans) );
TRANS_CHECK( trans_set_resol_lam(&trans, 20, 18, 2500.0, 2500.0) );
TRANS_CHECK( trans_set_trunc_lam(&trans, 9, 8) );
TRANS_CHECK( trans_setup(&trans) );
TRANS_CHECK( trans_inquire(&trans,"nvalue,mvalue") );
// Number of fields
int nscalar = 2;
int nvordiv = 1;
int nfld = 2*nvordiv+nscalar;
// Allocate test data
double* rgp1 = calloc( nfld * trans.ngptot, sizeof(double) );
double* rspvor1 = calloc( nvordiv * trans.nspec2, sizeof(double) );
double* rspdiv1 = calloc( nvordiv * trans.nspec2, sizeof(double) );
double* rspscalar1 = calloc( nscalar * trans.nspec2, sizeof(double) );
double* rmeanu1 = calloc( nvordiv, sizeof(double) );
double* rmeanv1 = calloc( nvordiv, sizeof(double) );
double* rgp2 = calloc( nfld * trans.ngptot, sizeof(double) );
double* rspvor2 = calloc( nvordiv * trans.nspec2, sizeof(double) );
double* rspdiv2 = calloc( nvordiv * trans.nspec2, sizeof(double) );
double* rspscalar2 = calloc( nscalar * trans.nspec2, sizeof(double) );
double* rmeanu2 = calloc( nvordiv, sizeof(double) );
double* rmeanv2 = calloc( nvordiv, sizeof(double) );
// Create random grid-point fields
for(int j=0; j<nfld; ++j)
for( int i=0; i<trans.ngptot; ++i )
rgp1[j*trans.ngptot+i] = randomDouble();
// Create random spectral fields
for( int i=0; i<trans.nspec2; ++i ) {
rspvor2[i*nvordiv+0] = randomDouble(); // vorticity
rspdiv2[i*nvordiv+0] = randomDouble(); // divergence
rspscalar2[i*nscalar+0] = randomDouble(); // scalar 1
rspscalar2[i*nscalar+1] = randomDouble(); // scalar 2
}
// Create random mean wind
for( int j=0; j<nvordiv; ++j) {
rmeanu2[j] = randomDouble();
rmeanv2[j] = randomDouble();
}
// Apply direct transform
struct DirTrans_t dirtrans = new_dirtrans(&trans);
dirtrans.nscalar = nscalar;
dirtrans.nvordiv = nvordiv;
dirtrans.rgp = rgp1;
dirtrans.rspscalar = rspscalar1;
dirtrans.rspvor = rspvor1;
dirtrans.rspdiv = rspdiv1;
dirtrans.rmeanu = rmeanu1;
dirtrans.rmeanv = rmeanv1;
TRANS_CHECK( trans_dirtrans(&dirtrans) );
// Compute spectral dot product
double adj_value_sp = 0.0;
for( int i=0; i<trans.nspec2; ++i ) {
adj_value_sp += rspvor1[i] * rspvor2[i];
adj_value_sp += rspdiv1[i] * rspdiv2[i];
}
for( int j=0; j<nscalar; ++j)
for( int i=0; i<trans.nspec2; ++i ) {
adj_value_sp += rspscalar1[i*nscalar+j] * rspscalar2[i*nscalar+j];
}
for( int j=0; j<nvordiv; ++j) {
adj_value_sp += rmeanu1[j] * rmeanu2[j];
adj_value_sp += rmeanv1[j] * rmeanv2[j];
}
printf("Spectral product = %.12f\n", adj_value_sp);
// Apply adjoint of the direct transform
struct DirTransAdj_t dirtrans_adj = new_dirtrans_adj(&trans);
dirtrans_adj.nscalar = nscalar;
dirtrans_adj.nvordiv = nvordiv;
dirtrans_adj.rspscalar = rspscalar2;
dirtrans_adj.rspvor = rspvor2;
dirtrans_adj.rspdiv = rspdiv2;
dirtrans_adj.rmeanu = rmeanu2;
dirtrans_adj.rmeanv = rmeanv2;
dirtrans_adj.rgp = rgp2;
TRANS_CHECK( trans_dirtrans_adj(&dirtrans_adj) );
// Compute grid-point dot product
double adj_value_gp = 0.0;
for(int j=0; j<nfld; ++j)
for( int i=0; i<trans.ngptot; ++i )
adj_value_gp += rgp1[j*trans.ngptot+i] * rgp2[j*trans.ngptot+i];
printf("Grid-point product = %.12f\n", adj_value_gp);
// Compare dot products
ASSERT( fabs(adj_value_sp - adj_value_gp ) / fabs( adj_value_gp ) < adjoint_tol );
// Deallocate arrays
free(rgp1);
free(rspvor1);
free(rspdiv1);
free(rspscalar1);
free(rmeanu1);
free(rmeanv1);
free(rgp2);
free(rspvor2);
free(rspdiv2);
free(rspscalar2);
free(rmeanu2);
free(rmeanv2);
TRANS_CHECK( trans_delete(&trans) );
}
// ----------------------------------------------------------------------------
int main ( int arc, char **argv ) {
trans_use_mpi( test_use_mpi() );
setbuf(stdout,NULL); // unbuffered stdout
// The adjoint test works for standard gaussian latitude grid
// with no points on the equator or poles.
// nsmax = nlat - 1
printf("-----------------------------\n");
test_lam_dirtrans_adjoint();
TRANS_CHECK( trans_finalize() );
return 0;
}
|
