-
Notifications
You must be signed in to change notification settings - Fork 0
/
primitives.c
executable file
·115 lines (72 loc) · 2.07 KB
/
primitives.c
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
#include <iostream>
#include <cmath>
#include <vector>
#include <stdio.h>
#include <functional>
#include "init.c"
double F(double Lor, double q2, double e, double d){
return q2/(1.-1./Lor/Lor)-pow((e*gamma_par*Lor*Lor-d*Lor)/(gamma_par*Lor*Lor-1.0),2);
}
double Fd(double Lor, double q2, double e, double d){
return -2*q2/pow((1.-1./Lor/Lor),2)*pow(Lor,-3)-2*((e*gamma_par*Lor*Lor-d*Lor)/(gamma_par*Lor*Lor-1.0))*((2*e*gamma_par*Lor-d)*(gamma_par*Lor*Lor-1.0)-(e*gamma_par*Lor*Lor-d*Lor)*(2*gamma_par*Lor))/pow((gamma_par*Lor*Lor-1.0),2);
}
// Quantities for Solution (8,9,10 delzana)
double W,E,D;
double Q2;
//error for the lorentz factor
double error = 0.01;
// how many itteration for N-R numerical solution
int n = 1;
int MAXITER = 100;
double x_lor;
int p;
// conversion from conservatives to primitves needed for each step
int conv_to_prims(ndvector main_grid){
for (int i = 2; i < main_grid.i; ++i){
for (int j = 0; j < main_grid.j; ++j)
{
for (int k = 0; k < main_grid.k; ++k)
{
p = i*main_grid.j*main_grid.k+j*main_grid.k+k; // fix the pointer to the "array"
// Build W,D,Q,E from (8,9,10) Del_Zana
E = main_grid.conservables_[p].e_;
D = main_grid.conservables_[p].rho_;
Q2=0;
for (int ii = 0; ii < n_comp; ++ii)
{
Q2 = Q2+pow(main_grid.conservables_[p].p_[ii],2);
}
n = 1;
x_lor = 1.2;
while( ( fabs(F(x_lor, Q2, E, D)) > error ) && ( n <= MAXITER ) )
{
x_lor = x_lor - ( F(x_lor, Q2, E, D)/ Fd(x_lor, Q2, E, D) );
n++;
}
std::cout << x_lor << ',' <<'\n';
}
}
}
return 0;
}
int main()
{
startup();
conv_to_prims(main_grid);
for (int i = 0; i < idim; ++i)
{
//std::cout << main_grid.array[i] << ' ' << '\n';
}
return 0;
}
/*
// x,y,z from grid.h used as reference to the real numerical distance rather than cell ID
main_grid.array[p].p_th = x[i];
main_grid.array[p].rho = y[j];
for (int ii = 0; ii < ndim; ++ii)
{
main_grid.array[p].u[ii] = x[i]/(x[i]+z[k]);
}
// builed conservable
main_grid.conservables_[p].build(main_grid.array[p]);
*/