Browse files

Using Fortran order in NumPy, changed tabs to spaces, other fixes

  • Loading branch information...
1 parent 8c2eb46 commit 9d09270ef2b1e6ca5c1cd7c73bdb84fef7f00d59 @mbarriault committed Nov 17, 2011
Showing with 139 additions and 123 deletions.
  1. +22 −19 Pulsar/main.cpp
  2. +3 −3 T3nsors/Partial.cpp
  3. +3 −2 T3nsors/Stream.cpp
  4. +1 −1 T3nsors/Stream.h
  5. +1 −1 T3nsors/System.cpp
  6. +2 −1 T3nsors/System.h
  7. +7 −4 bin/t3.l2plot.py
  8. +34 −28 bin/t3.quadfit.py
  9. +66 −64 t3py/__init__.py
View
41 Pulsar/main.cpp
@@ -29,11 +29,11 @@ Vector& dB2 = (Vector&)dx2.at(2);
#include <cmath>
using namespace T3;
-Params PulsarParams(real, real, real, real, real, real);
+Params PulsarParams(real, real, real, real, real, real, real);
class Pulsar : public System {
public:
- Axisymmetric Del;
+ Axisymmetric Nabla;
Partial& r;
Partial& theta;
Params P;
@@ -45,14 +45,15 @@ class Pulsar : public System {
Stream EB, divB;
Stream S;
- Pulsar(int n, real s, int rot, std::string iid, Params P) : P(P), Del(Axisymmetric(n, 1., P["RLC"]/P["Omega"])), r(Del[0]) , theta(Del[1]) {
+ Pulsar(int n, real s, int rot, std::string iid, Params P) : P(P), Nabla(Axisymmetric(n, P["Rs"]/P["Omega"], P["RLC"]/P["Omega"])), r(Nabla[0]) , theta(Nabla[1]) {
+ Del = &Nabla;
real T = rot*2*M_PI/P["Omega"];
- System::t = Partial(0, 0., s*Del[0].d, T, this);
- Stream::Del = &Del;
+ t = Partial(0, 0., s*Nabla[0].d, T, this);
+ //Stream::Del = &Del;
id = mkdir(iid + timecoord());
- epsilon = 0.5;
+ epsilon = 0.5/Nabla[0].d;
c = P["alpha"]/(t(-1) * P["mu"]*P["mu"]/pow(r(0),6.) * (1 + pow(P["Omega"]*r(0), 2.) ) );
- push_back(Stream(Scalar(), "u", this));
+ push_back(Stream(Scalar(), "Psi", this));
push_back(Stream(Vector(), "E", this));
push_back(Stream(Vector(), "B", this));
@@ -67,6 +68,7 @@ class Pulsar : public System {
void Initialize() { SLINK
Set x = Source(*this, Tensor::N[0]);
+ Psi = x[0];
E = x[1];
B = x[2];
@@ -81,6 +83,7 @@ class Pulsar : public System {
Set Source(Set x, int in=1) { LINK
PFRO(i,0,in) FOR(j,Tensor::N[1]) FOR(k,Tensor::N[2]) {
+ Psi(i,j,k) = 0.;
B(0)(i,j,k) = 2*P["mu"]*cos(theta(j)) / pow(r(i),3.);
B(1)(i,j,k) = P["mu"]*sin(theta(j)) / pow(r(i),3.);
B(2)(i,j,k) = 0.;
@@ -102,7 +105,7 @@ class Pulsar : public System {
}
Scalar CalcRho(Set x) { LINK
- Scalar rho = Del*E;
+ Scalar rho = Nabla*E;
FOR(o,Tensor::N.Pr()) rho[o] /= 4*M_PI;
return rho;
}
@@ -127,17 +130,17 @@ class Pulsar : public System {
}
Scalar CalcDivB(Set x) { LINK
- return Del*B;
+ return Nabla*B;
}
Set RHS(Set& x) { LINK
x = Source(x);
Set dx = x; DLINK
Scalar divB = CalcDivB(x);
Vector J = CalcJ(x);
- Vector curlB = (Del&B);
- Vector curlE = (Del&E);
- Vector gradPsi = Del(Psi);
+ Vector curlB = (Nabla&B);
+ Vector curlE = (Nabla&E);
+ Vector gradPsi = Nabla(Psi);
Vector F1;
Vector F2;
@@ -154,14 +157,14 @@ class Pulsar : public System {
dB = LC(-1., &curlE, 1., &gradPsi, -c, &F2, NULL);
Set dx2 = dx; D2LINK
- PFOR(j,Tensor::N[1]) FOR(k,Tensor::N[2]) {
+ /*PFOR(j,Tensor::N[1]) FOR(k,Tensor::N[2]) {
dPsi2(0,j,k) = 0.5*( dPsi(0,j,k) + dB(0)(0,j,k) );
dE2(1)(0,j,k) = 0.5*( dE(1)(0,j,k) - dB(2)(0,j,k) );
dE2(2)(0,j,k) = 0.5*( dE(2)(0,j,k) + dB(1)(0,j,k) );
dB2(0)(0,j,k) = 0.5*( dB(0)(0,j,k) + dPsi(0,j,k) );
dB2(1)(0,j,k) = 0.5*( dB(1)(0,j,k) + dE(2)(0,j,k) );
dB2(2)(0,j,k) = 0.5*( dB(2)(0,j,k) - dE(1)(0,j,k) );
- }
+ }*/
PFOR(j,Tensor::N[1]) FOR(k,Tensor::N[2]) {
dPsi2(-1,j,k) = 0.5*( dPsi(-1,j,k) - dB(0)(-1,j,k) );
dE2(1)(-1,j,k) = 0.5*( dE(1)(-1,j,k) + dB(2)(-1,j,k) );
@@ -209,8 +212,9 @@ class Pulsar : public System {
}
};
-Params PulsarParams(real RLC, real k, real mu, real Omega, real alpha, real zeta) {
+Params PulsarParams(real Rs, real RLC, real k, real mu, real Omega, real alpha, real zeta) {
Params P;
+ P["Rs"] = Rs;
P["RLC"] = RLC;
P["k"] = k;
P["mu"] = mu;
@@ -222,10 +226,9 @@ Params PulsarParams(real RLC, real k, real mu, real Omega, real alpha, real zeta
int main (int argc, const char * argv[])
{
- Pulsar(48, 0.125, 2, "pulsar-test", PulsarParams(25., 1., 1., 0.75, 0., 1.)).Run();
-/* for ( real Omega = 0.5; Omega < 0.61; Omega += 0.01 )
- Pulsar(12, 0.125, 1, "pulsar-test", PulsarParams(1., 1., Omega, 0., 1.)).Run();*/
- // insert code here...
+ //Pulsar(18, 0.25, 5, "pulsar-test", PulsarParams(0.2, 4., 1., 1., 0.5, 100., 0.)).Run();
+ for ( real Omega = 0.5; Omega < 0.61; Omega += 0.01 )
+ Pulsar(18, 0.25, 1, "pulsar-static", PulsarParams(0.2, 4., 1., 1., Omega, 0., 1.)).Run();
return 0;
}
View
6 T3nsors/Partial.cpp
@@ -44,9 +44,9 @@ T3::Partial T3::Partial::Cartesian(int p, real d, real a, real b, Object* parent
T3::Partial T3::Partial::Azimuth(int p, int n, Object* parent) {
Partial D(p, 0., M_PI, n, parent);
- D.a += D.d;
- D.b -= D.d;
- D.n -= 1;
+ D.a += D.d/2;
+ D.b -= D.d/2;
+// D.n -= 1;
D.bnd = bounds_bounce;
return D;
}
View
5 T3nsors/Stream.cpp
@@ -11,8 +11,8 @@
#include <iostream>
int T3::Stream::nt_min(3);
-int T3::Stream::nt_max(5);
-T3::Connection* T3::Stream::Del(NULL);
+int T3::Stream::nt_max(20);
+//T3::Connection* T3::Stream::Del(NULL);
T3::Stream::Stream() {
@@ -21,6 +21,7 @@ T3::Stream::Stream() {
T3::Stream::Stream(T3::Tensor x, std::string id, Object* parent) {
this->parent = parent;
this->id = id;
+ Del = ((System*)parent)->Del;
x.t = 0;
x.parent = this;
push_back(x);
View
2 T3nsors/Stream.h
@@ -21,7 +21,7 @@ namespace T3 {
public:
static int nt_min;
static int nt_max;
- static Connection* Del;
+ Connection* Del;
H5::H5File file;
H5::Group datagrp;
View
2 T3nsors/System.cpp
@@ -8,7 +8,7 @@
#include "System.h"
-T3::Partial T3::System::t(Partial::Cartesian(0, 10, 0., 1.));
+//T3::Partial T3::System::t(Partial::Cartesian(0, 10, 0., 1.));
void T3::System::Evolve() {
PreEvolve();
View
3 T3nsors/System.h
@@ -18,7 +18,8 @@ namespace T3 {
void Evolve();
Set Dissipate(const Set&);
public:
- static Partial t;
+ Partial t;
+ Connection* Del;
real epsilon;
virtual Set RHS(Set&) = 0;
virtual int Condition() = 0;
View
11 bin/t3.l2plot.py
@@ -11,10 +11,13 @@
l2f = []
theta,r = numpy.meshgrid(F.x[1],F.x[0])
for f in F:
- val = 0
- f = f.reshape(r.shape)
- g = f*(r**2)*numpy.sin(theta)
- l2f.append( math.sqrt(numpy.dot(g.flat,g.flat)/g.size) )
+ val = 0
+ try:
+ f = f.reshape(r.shape)
+ except ValueError:
+ f = f[0].reshape(r.shape)
+ g = f*(r**2)*numpy.sin(theta)
+ l2f.append( math.sqrt(numpy.dot(g.flat,g.flat)/g.size) )
pyplot.plot(F.t, l2f)
print l2f
View
62 bin/t3.quadfit.py
@@ -10,34 +10,40 @@
val = sys.argv[2]
runs = sys.argv[3:]
-Ms = numpy.linspace(0.8, 1., 3)
+Ms = numpy.linspace(0.2, 1., 9)
+#Ms = (0.5,1.)
for M in Ms:
- V = []
- for run in runs:
- F = t3py.hdf(os.path.join(run,func))
- v = F.getparam(val,os.path.join(run,"info.txt"))
- S = F.shell(-1)
- R = F.x[0][0]
- n = F.x[0][-1]*v
- m = M*n
- i = int( (m - R*v)/(n-R*v) * (len(F.x[0]) - 1) )
- V.append( (v, S[i]) )
-
- v = [x[0] for x in V]
- L = [x[1] for x in V]
- Lfit = numpy.polyfit(v, L, 4)
-
- h = (M - Ms[0])/(Ms[-1]-Ms[0]) * 2./3
- c = colors.hsv_to_rgb(numpy.array([[[h,1,1]]]))
- print Lfit
-
- C = numpy.zeros((len(L),3))
- for i,Cc in enumerate(C):
- C[i] = c[0,0]
-
- pyplot.xlabel(val)
- pyplot.ylabel(func[:-4])
- pyplot.scatter(v, L, color=C)
- pyplot.plot(v, numpy.poly1d(Lfit)(v), color=c[0,0], label=str(M))
+ V = []
+ for run in runs:
+ F = t3py.hdf(os.path.join(run,func))
+ v = F.getparam(val,os.path.join(run,"info.txt"))
+ Sr = numpy.zeros(F.x[0].shape)
+ for i,t in enumerate(F.t):
+ Sr += F.shell(i,0)
+ Sr /= len(F.t)
+ #Sr = F.shell(-1,0)
+ R = F.x[0][0]
+ n = F.x[0][-1]*v
+ m = M*n
+ i = int( (m - R*v)/(n-R*v) * (len(F.x[0]) - 1) )
+ V.append( (v, Sr[i]) )
+
+ v = [x[0] for x in V]
+ L = [x[1] for x in V]
+ Lfit = numpy.polyfit(v, L, 4)
+
+ h = (M - Ms[0])/(Ms[-1]-Ms[0]) * 2./3
+ c = colors.hsv_to_rgb(numpy.array([[[h,1,1]]]))
+ print Lfit
+ print L
+
+ C = numpy.zeros((len(L),3))
+ for i,Cc in enumerate(C):
+ C[i] = c[0,0]
+
+ pyplot.xlabel(val)
+ pyplot.ylabel(func[:-4])
+ pyplot.scatter(v, L, color=C)
+ pyplot.plot(v, numpy.poly1d(Lfit)(v), color=c[0,0], label=str(M))
pyplot.legend()
pyplot.show()
View
130 t3py/__init__.py
@@ -3,67 +3,69 @@
import math
class hdf:
- def __init__(self, file):
- self.file = file
- self.f = h5py.File(file, 'r')
- self.t = sorted(self.f['data'].keys(), key=lambda x : float(x))
- self.cnt = 0
- F = self.f['data'][self.t[0]]
- if type(F) is h5py._hl.dataset.Dataset:
- self.rank = 0
- self.shape = F.shape
- elif type(F) is h5py._hl.group.Group:
- self.rank = len(F.keys())
- self.shape = F['0'].shape
- self.x = []
- X = sorted(self.f['coords'].items())
- for i,(c,x) in enumerate(X):
- ar = numpy.arange(x[0], x[1]+x[2]/2, x[2])
- if len(ar) > self.shape[i]:
- ar = ar[:-1]
- self.x.append(ar)
- def __del__(self):
- self.f.close()
- def __getitem__(self, index):
- key = self.t[index]
- F = self.f['data'][key]
- if self.rank == 0:
- f = numpy.array(F)
- else:
- f = []
- for c in self.f['data'][key]:
- f.append(numpy.array(F[c]))
- f = numpy.array(f)
- return f
- def size(self):
- return len(self.t)
- def __iter__(self):
- return self
- def next(self):
- try:
- result = self[self.cnt]
- self.cnt += 1
- return result
- except:
- raise StopIteration
- def shell(self,i,c=0):
- if self.rank == 0:
- f = self[i]
- else:
- f = self[i][c]
- shape = f.shape[0:1]
- s = numpy.zeros(shape)
- for j in xrange(f.shape[1]):
- for k in xrange(f.shape[2]):
- s += f[:,j,k] * (self.x[0])**2 * numpy.sin(self.x[1][j])
- s -= 0.5 * ( f[:,0,0] * (self.x[0])**2 * numpy.sin(self.x[1][0]) + f[:,-1,0] * (self.x[0])**2 * numpy.sin(self.x[1][-1]) )
- s *= 2*math.pi * (self.x[1][-1]-self.x[1][0])/self.x[1].size
- return s
- def getparam(self,param,file="info.txt"):
- info = open(file, 'r')
- for line in info.readlines():
- if param in line:
- val = float( line.split(' ')[2] )
- break
- return val
-
+ def __init__(self, file):
+ self.file = file
+ self.f = h5py.File(file, 'r')
+ self.t = sorted(self.f['data'].keys(), key=lambda x : float(x))
+ self.cnt = 0
+ F = self.f['data'][self.t[0]]
+ if type(F) is h5py._hl.dataset.Dataset:
+ self.rank = 0
+ self.shape = F.shape
+ elif type(F) is h5py._hl.group.Group:
+ self.rank = len(F.keys())
+ self.shape = F['0'].shape
+ self.x = []
+ X = sorted(self.f['coords'].items())
+ for i,(c,x) in enumerate(X):
+ ar = numpy.arange(x[0], x[1]+x[2]/2, x[2])
+ if len(ar) > self.shape[i]:
+ ar = ar[:-1]
+ self.x.append(ar)
+ def __del__(self):
+ self.f.close()
+ def __getitem__(self, index):
+ key = self.t[index]
+ F = self.f['data'][key]
+ if self.rank == 0:
+ f = numpy.array(F)
+ f = f.reshape(f.shape, order='F')
+ else:
+ f = []
+ for c in self.f['data'][key]:
+ Fc = numpy.array(F[c])
+ Fc = Fc.reshape(Fc.shape, order='F')
+ f.append(Fc)
+ return f
+ def size(self):
+ return len(self.t)
+ def __iter__(self):
+ return self
+ def next(self):
+ try:
+ result = self[self.cnt]
+ self.cnt += 1
+ return result
+ except:
+ raise StopIteration
+ def shell(self,i,c=0):
+ if self.rank == 0:
+ f = self[i]
+ else:
+ f = self[i][c]
+ shape = f.shape[0:1]
+ s = numpy.zeros(shape)
+ for j in xrange(f.shape[1]):
+ for k in xrange(f.shape[2]):
+ s += f[:,j,k] * (self.x[0])**2 * numpy.sin(self.x[1][j])
+ s -= 0.5 * ( f[:,0,0] * (self.x[0])**2 * numpy.sin(self.x[1][0]) + f[:,-1,0] * (self.x[0])**2 * numpy.sin(self.x[1][-1]) )
+ s *= 2*math.pi * (self.x[1][-1]-self.x[1][0])/self.x[1].size
+ return s
+ def getparam(self,param,file="info.txt"):
+ info = open(file, 'r')
+ for line in info.readlines():
+ if param in line:
+ val = float( line.split(' ')[2] )
+ break
+ return val
+

0 comments on commit 9d09270

Please sign in to comment.