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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
#include <iostream>
#include <math.h>
#include <fstream>
#include <algorithm>
using namespace std;
const double MM = 8;
const double Rin = 1.0;
const double Rout = 2.0;
const double Zmin = 0.0;
const double Zmax = 3.14159;
const double tend = 0.04;
double dtout = 1.0;
const double D = 0.1;
const double factor = 0.96;
double dr = 1.0/MM;
double dz = 1.0/MM; // note right now assuming R and Z mesh have same MM
const int MR = (Rout-Rin)*MM;
const int MZ = (Zmax-Zmin)*MM;
double V[MR+1][MZ+1] = {0.0};
double r[MR+1] = {0.0};
double z[MZ+1] = {0.0};
double Fr[MR+1][MZ+1] = {0.0};
double Fz[MR+1][MZ+1] ={0.0};
double Ar[MR+1][MZ+1] = {0.0};
double Az[MR+1][MZ+1] = {0.0};
double dtEXPL = 1.0/(2.0*D*(1.0/(dr*dr)+1.0/(dz*dz)));
double dt = factor*dtEXPL;
double MaxSteps = int((tend/dt)+1.0);
const double pi = 4.0*atan(1.0);
double Texact[MR+1][MZ+1] = {0.0};
double T[MR+1][MZ+1] = {0.0};
double time_var = 0.0;
double Err[MR+1]; ////what size?
//double sintest;
double EXACT(double Texact[MR+1][MZ+1], double r[], double z[], double time_var, int MR, int MZ) {
for(int i = 1 ; i < MR+1; ++i ){
for(int j = 1 ; j < MZ+1; ++j ){
Texact[i][j] = exp(-time_var)*log(r[i])*sin(z[j]);
}
}
}
double FLUX(double time_var, double r[], double z[], double Ar[][MZ+1], double Az[][MZ+1],
double T[][MZ+1], double Texact[][MZ+1], double Fr[][MZ+1], double Fz[][MZ+1],
double D, double dr, double dz, int MR, int MZ){
//boundary cnditions
for(int i = 0; i<MR+1; ++i) {
Texact[i][0] = 0;
Texact[i][MZ] = 0;
T[i][0] = 0;
T[i][MZ] = 0;
Fz[i][1] = -D*(T[i][1] - T[i][0])/(z[1]-z[0]);
Fz[i][MZ] = -D*(T[i][MZ] - T[i][MZ-1])/(z[MZ]-z[MZ-1]);
cout << Ar[2][2] << "flux1 " << endl;
}
for(int j = 0; j<MZ+1; ++j) {
Texact[0][j] = 0;
Texact[MR][j] = exp(-time_var)*log(2.0)*sin(z[j]);
T[0][j] = 0;
T[MR][j] = exp(-time_var)*log(2.0)*sin(z[j]);
Fz[1][j] = -D*(T[1][j] - T[0][j])/(r[1]-r[0]);
Fz[MR][j] = -D*(T[MR][j] - T[MR-1][j])/(r[MR]-r[MR-1]);
}
// flux calc
for(int i = 2; i<MZ+1; ++i){
for(int j = 2; j<MZ+1; ++j){
if (j < 9) {
cout << Ar[0][0] << " " << Ar[0][1] << " flux32 " << i << " i " << j <<" j "<< endl;}
Fr[i][j] = -D*(T[i][j] - T[i-1][j])/dr;
if (j < 9) {
cout << Ar[0][0] << " " << Ar[0][1] << " flux33 " << i << " i " << j <<" j "<< endl;}
Fz[i][j] = -D*(T[i][j] - T[i][j-1])/dz;
}
}
}
double PDE(double T[][MZ+1], double Fr[][MZ+1], double Fz[][MZ+1], double dt,
double time_var, double Ar[][MZ+1], double Az[][MZ+1], double V[][MZ+1], int MR, int MZ){
for(int i = 1; i<MR+1; ++i){
cout << Ar[2][2] << "pde1 "<< endl;
for(int j = 1; j<MZ+1; ++j){
T[i][j] = T[i][j] + (dt/V[i][j])*((Fr[i][j]*Ar[i][j])-(Fr[i+1][j]*Ar[i+1][j])+(Fz[i][j]*Az[i][j])-(Fz[i][j+1]*Az[i][j+1])) ;
}
}
}
double OUTPUT(double time_var, double r[], double z[], double T[][MZ+1],double Fr[][MZ+1],
double Fz[][MZ+1], double dr, double dz, double Texact[][MZ+1], double Ar[][MZ+1],
double Az[][MZ+1], double Err[], int MZ, int MR) {
ofstream myfile;
myfile.open ("lab5out.txt",ios::app); // opens in mode to append
myfile << "dr, dz " << dr << " " << dz <<"\n" ; // headers
myfile << "#time r z T Texact FR FZ Ar Az: " <<"\n" ; // headers
for(int m=0; m<MR+2; ++m){
for(int n=0; n<MZ+2; ++n){
myfile << time_var << " " << r[m] << " " << z[n] << " "<< T[m][n]
<< " "<< Texact[m][n] <<" "<< Fr[m][n] <<" " <<Fz[m][n] << " "<<
Ar[m][n] << " " << Az[m][n] << " m= " << m << " n= " << n << "\n"; // report all values for current time - report time, x position and U
}}
myfile <<"\n";// skip a line after each time
myfile.close();
}
int main() {
//init mesh for r
r[0] = Rin;
r[1] = Rin + dr/2;
r[MR+1] = Rout;
for (int i=2; i<MR+1 ; ++i) {
r[i] = r[1] + (i-1)*dr;
}
//init mesh for z
z[0] = Zmin;
z[1] = Zmin + dz/2;
z[MZ+1] = Zmax;
for (int i=2; i<MZ+1 ; ++i) {
z[i] = z[1] + (i-1)*dz;
}
//build geometry for areas and volume
for(int i = 0; i<MR+1; ++i) {
for(int j = 0; j<MZ+1; ++j){
Ar[i][j] = 2.0*pi*r[i]*dz;
Az[i][j] = 2.0*pi*r[i]*dr;
}
}
for(int i = 0; i<MR+1; ++i) {
for(int j = 0; j<MZ+1; ++j){
V[i][j] = 2.0*pi*r[i]*dr*dz;
}
}
// initial condition
for(int i = 1; i<MR+1; ++i) {
for(int j = 1; j<MZ+1; ++j){
Texact[i][j] = log(r[i])*sin(z[j]);
T[i][j] = log(r[i])*sin(z[j]);
}
}
ofstream myfile; // open file to clear it for later use, and plot the initial
myfile.open ("lab5out.txt",ios::trunc);
myfile << "dr, dz, dt " << dr << " " << dz << " " << dt <<"\n" ; // headers
myfile << "#time r z T Texact FR FZ Ar Az: " <<"\n" ; // headers
for(int m=0; m<MR+1; ++m){
for(int n=0; n<MZ+1; ++n){
myfile << time_var << " " << r[m] << " " << z[n] << " "<< T[m][n] << " "<< Texact[m][n] <<" "<< Fr[m][n] <<" " <<Fz[m][n] << " "<< Ar[m][n] << " " << Az[m][n] << " m= " << m << " n= " << n << "\n"; // report all values for current time - report time, x position and U
}
}
myfile <<"\n";// skip a line after each time
myfile.close();
//// inside time step
int nsteps = 0;
dtout = max(dtout,dt);
int tout = dtout;
// time stepping
for (int n = 1; n < MaxSteps; n++) {
time_var = n*dt;
FLUX(time_var,r,z, Ar, Az, T, Texact, Fr, Fz, D, dr, dz, MR, MZ);
PDE(T,Fr,Fz,dt,time_var,Ar,Az, V, MR, MZ);
EXACT(Texact, r, z, time_var, MR, MZ);
if (time_var >= dtout){
OUTPUT(time_var,r,z,T,Fr,Fz,dr,dz,Texact,Ar,Az,Err,MZ,MR);
dtout = tout + dtout;
}
}
}
|
