Skip to content

Commit

Permalink
Dodano funkcje do raportowania nieprawidłowych wartości liczbowych.
Browse files Browse the repository at this point in the history
  • Loading branch information
Peter committed Nov 20, 2011
1 parent c671e9a commit 01d3b49
Show file tree
Hide file tree
Showing 4 changed files with 139 additions and 45 deletions.
5 changes: 4 additions & 1 deletion .gitignore
View file
@@ -1,4 +1,7 @@
#Ignoruj pliki obiektów
#Ignoruj pliki obiektów, tymczasowe pliki edytora i skrypty octave
*.h
*.gch
*.o
*.oct
*.m
*.swp
10 changes: 5 additions & 5 deletions src/NLS_solver.cpp
View file
Expand Up @@ -82,17 +82,17 @@ DEFUN_DLD( NLS_solver, args, nargout , "Non-linear Schrodinger equation solver."
//punkt czasowy z którego zaczynamy


complex x;
std::vector<complex> in_psi0( 0);
std::vector<complex> in_V(0);
std_complex x;
std::vector<std_complex> in_psi0( 0);
std::vector<std_complex> in_V(0);
std::vector<double> in_n_r0(0);
std::vector<double> in_P_l(0);
physical_constants in_consts;

for(int i =0; i< spatial_size ;i++)
{
in_psi0.push_back( complex( psi0(i).real(), psi0(i).imag() ));
in_V.push_back( complex( V(i).real(), V(i).imag() ));
in_psi0.push_back( std_complex( psi0(i).real(), psi0(i).imag() ));
in_V.push_back( std_complex( V(i).real(), V(i).imag() ));
in_n_r0.push_back( n_r0(i) );
in_P_l.push_back( P_l(i) );
}
Expand Down
108 changes: 83 additions & 25 deletions src/solve.cpp
View file
Expand Up @@ -2,7 +2,7 @@

const double PI = 3.1415926535;

solution::solution(const std::vector<complex>& temp, const std::vector<complex>& temp_V,
solution::solution(const std::vector<std_complex>& temp, const std::vector<std_complex>& temp_V,
const std::vector<double>& temp_n, const std::vector<double>& temp_p,
double xstep =0.01, double tstep=0.001, physical_constants C = physical_constants() )
{
Expand All @@ -19,20 +19,25 @@ solution::solution(const std::vector<complex>& temp, const std::vector<complex>&
CONST.R = C.R;

data_size = temp.size();
V = std::vector<complex>(temp_V);
V = std::vector<std_complex>(temp_V);
n_r = std::vector<double>(temp_n);
P_l = std::vector<double>(temp_p);

//inicjacja fftw
FT=fourier_transform();
psi = reinterpret_cast<complex*>( fftw_malloc(sizeof(complex) * data_size) );

//FIXME -- czy ta konwersja jest bezpieczna? znaleźć inny sposób
psi = reinterpret_cast<std_complex*>( fftw_malloc(sizeof(std_complex) * data_size) );
FT.init( psi, psi, data_size);

for(int i = 0; i < data_size; ++i)
{
psi[i] = temp[i];
}

error_state = 0;
step = 0;

}//spoko


Expand All @@ -44,11 +49,14 @@ solution::solution(const std::vector<complex>& temp, const std::vector<complex>&
//
// 2) dN_dt = Pl (x,t ) − γR n_R (x,t ) − R(n_R (x,t ))|ψ(x,t )|2
//


solution::~solution()
{
fftw_free(psi);
//fftw_free(A_n);
}
//spoko loko


void solution::first_step()
{
Expand All @@ -58,7 +66,7 @@ void solution::first_step()
for(int i =0; i<data_size ; ++i)
{
//normalizacja
psi[i] = psi[i] / complex(data_size,0.0);
psi[i] = psi[i] / std_complex(data_size,0.0);
//coś co wyskakuje przy liczeniu pochodnej
//numer indeksu po którym częstotliwość zaczyna być dodatnia
int i_max = (data_size-data_size%2)/2 - 1 + data_size%2;
Expand All @@ -69,7 +77,7 @@ void solution::first_step()
}
double k_square = std::pow( PI*k*2.0 / (data_size*step_spatial), 2.0);
//std::cout<<std::sqrt(k_square)<<std::endl;
psi[i] = psi[i] * std::exp(complex(0.0,-k_square ) * step_temporal * 0.5*CONST.h_bar/CONST.m ) ;
psi[i] = psi[i] * std::exp(std_complex(0.0,-k_square ) * step_temporal * 0.5*CONST.h_bar/CONST.m ) ;
}


Expand All @@ -78,18 +86,18 @@ void solution::first_step()

void solution::second_step()
{
complex phase(0.0,0.0);
std_complex phase(0.0,0.0);
for(int i = 0 ; i < data_size; ++i)
{
//EWOLUCJA "psi"
//ewolucja napędzane nieliniowością i potencjałem
phase = complex(0.0, CONST.g_c * std::norm(psi[i]) );
phase += complex(0.0,1.0) * V[i] ;
phase = std_complex(0.0, CONST.g_c * std::norm(psi[i]) );
phase += std_complex(0.0,1.0) * V[i] ;

//ewolucja napędzana przez rozpraszanie i zmiany w ilości polarytonów
phase += complex( 0.5* (CONST.lam_c - R(n_r[i]))*CONST.h_bar, CONST.g_r *n_r[i] ); // plus działanie funkcji R na n_r
phase += std_complex( 0.5* (CONST.lam_c - R(n_r[i]))*CONST.h_bar, CONST.g_r *n_r[i] ); // plus działanie funkcji R na n_r
//aplikacja do funkcji falowej
psi[i] = psi[i] * std::exp( -phase * step_temporal /CONST.h_bar );
psi[i] = psi[i] * std::exp( -phase * (step_temporal /CONST.h_bar) );

//EWOLUCJA "n_r"
n_r[i] += ( P_l[i] - CONST.lam_r * n_r[i] - R(n_r[i] ) * std::norm(psi[i]) )*step_temporal;
Expand All @@ -104,19 +112,19 @@ void solution::make_time_step()
//jak wskazywałaby nazwa krok drugi jest wykonywany pierwszy
first_step();
second_step();
++step;
}//spoko

void solution::evolution(int steps )
{
for(int i = 0; i < steps; i++)
{
make_time_step();
apply_boundary_condition(complex(0.0,0.0), complex(0.0,0.0), 0.1 );
apply_boundary_condition(std_complex(0.0,0.0), std_complex(0.0,0.0), 0.1 );
}
}//spoko

//
void solution::apply_boundary_condition(complex a=complex(0.0,0.0), complex b=complex(0.0,0.0), double c = 0.1)
void solution::apply_boundary_condition(std_complex a=std_complex(0.0,0.0), std_complex b=std_complex(0.0,0.0), double c = 0.1)
{
int falloff = std::floor( c * static_cast<double>(data_size) );
int k;
Expand Down Expand Up @@ -150,31 +158,81 @@ double solution::abs_val( )
wynik+= std::norm(psi[i]) * step_spatial;
return std::sqrt(wynik);
}//spoko loko
const std::vector<complex> solution::output_psi()

const std::vector<std_complex> solution::output_psi()
{
std::vector<complex> out = std::vector<complex>(0);
std::vector<std_complex> out = std::vector<std_complex>(0);
for(int i = 0; i < data_size; ++i )
out.push_back(psi[i]);
return out;
}//spoko loko

const complex solution::output_psi(int i)
const std_complex solution::output_psi(int i)
{
return psi[i];
}
//spoko loko
}//spoko loko

const std::vector<double> solution::output_n_r()
{
std::vector<double> out = std::vector<double>(0);
for(int i = 0; i < data_size; ++i )
out.push_back(n_r[i]);
return out;
return n_r;
}//spoko loko

const double solution::output_n_r(int i)
{
return n_r[i];
}
}//spoko loko

const std::string solution::report_exceptions()
{
//FIXME -- czy takie operacje dają wynik niezależny od implementacji "int"
error_state = 0;
std::ostringstream report;

std::string comment("NAN value found in: ");

for(int i = 0;i<data_size; ++i)
{
if( nan_exception::is_nan(n_r[i]) && (error_state & nan_exception::nan_n_r) )
{
error_state = error_state | nan_exception::nan_n_r;
report << comment;
report << std::string("n_r");
report << std::endl;
}
if( nan_exception::is_nan(psi[i]) && (error_state & nan_exception::nan_psi) )
{
error_state = error_state | nan_exception::nan_psi;
report << comment;
report << std::string("psi");
report << std::endl;
}
if( nan_exception::is_nan(V[i]) && ( error_state & nan_exception::nan_V) )
{
error_state = error_state | nan_exception::nan_V;
report << comment;
report << std::string("V");
report << std::endl;
}
if( nan_exception::is_nan(P_l[i]) && (error_state & nan_exception::nan_P_l) )
{
error_state = error_state | nan_exception::nan_P_l;
report << comment;
report << std::string("P_l");
report << std::endl;
}
}

if(error_state == 0)
{
report << std::string("No exception::isnans found in step");
}

report << std::string("in step: ");
report << step;

return report.str();
}//spoko loko

//-------------------------PHYSICAL_CONSTANTS------------------------------
physical_constants::physical_constants()
{
Expand All @@ -184,7 +242,7 @@ physical_constants::physical_constants()
m = 0.5;
lam_c = 0.0;
lam_r = 0.0;
R=0.0;
R=0.0;
}

//-------------------FFTW3-------------------------------------------------
Expand Down
61 changes: 47 additions & 14 deletions src/solve.h
View file
Expand Up @@ -8,11 +8,15 @@

#include<vector>
#include<complex>

//funkcje matematyczne
#include<cmath>
//transformata fouriera
#include <fftw3.h>
//raportowanie błędów
#include<sstream>

//PRZESTRZEŃ NAZW - chyba nie potrzebna jak na razie
typedef std::complex<double> std_complex;


enum direction
Expand All @@ -21,12 +25,27 @@ enum direction
backward = 1
};

namespace nan_exception
{

//PRZESTRZEŃ NAZW - chyba nie potrzebna jak na razie
typedef std::complex<double> complex;
enum Enum
{
nan_psi = 1,
nan_n_r = 2,
nan_P_l = 4,
nan_V = 8,
};

bool inline is_nan( std_complex x )
{
if( std::isnan( x.real() ) || std::isnan( x.imag() ) )
return true;
else
return false;
}

}

//TYP DO MACIERZY LICZB RZECZYWISTYCH
//typedef std::vector< std::vector<double> > matrix_double;

//TRANSFORMACJA FOURIERA
class fourier_transform
Expand All @@ -39,7 +58,8 @@ class fourier_transform
fftw_plan forward_plan;
fftw_plan backward_plan;
public:
void inline init( complex* in , complex* out, int size)
//FIXME -- czy coś takiego rzutowanie jest bezpieczne, może da się lepiej?
void inline init( std_complex* in , std_complex* out, int size)
{
forward_plan = fftw_plan_dft_1d( size, reinterpret_cast<fftw_complex*>(in),
reinterpret_cast<fftw_complex*>(out), FFTW_FORWARD, FFTW_PATIENT);
Expand Down Expand Up @@ -75,24 +95,30 @@ class solution
private:
//zmienne przechowujące wartość rozwiązania dla ustalonego czasu "t"
//warunki brzegowe
//complex boundary_begin;
//complex boundary_end;
complex* psi;
//std_complex boundary_begin;
//std_complex boundary_end;

//FIXME - może jest bezpieczniejszy sposób żeby pogodzić fftw i wskaźniki?
std_complex* volatile psi;

std::vector<double> n_r;

//DANE RÓWNANIA
//długość kroku
double step_spatial;
double step_temporal; //krok czasowy może być urojony odpowiada to operatorowi nie
int data_size;
std::vector<complex> V;
std::vector<std_complex> V;
std::vector<double> P_l;

//STAŁE FIZYCZNE
physical_constants CONST;

fourier_transform FT;


//to zanczy dzielenie przez zero albo wyjście poza zakres

//FUNKCJA, która może być dostarczona z zewnątrz, ale na razie po prostu jest metodą
//żeby być bardziej "explicit"
double inline R( double s )
Expand All @@ -104,17 +130,24 @@ class solution
void first_step();
void second_step();
void make_time_step();
void apply_boundary_condition(complex start, complex end, double falloff);
void apply_boundary_condition(std_complex start, std_complex end, double falloff);

//DANE KONTROLNE
int error_state;
int step;

public:
solution(const std::vector<complex>& temp, const std::vector<complex>& temp_V,
solution(const std::vector<std_complex>& temp, const std::vector<std_complex>& temp_V,
const std::vector<double>& temp_n, const std::vector<double>& temp_p,
double xstep, double tstep, physical_constants C );
~solution();
void evolution(int steps );
const std::vector<complex> output_psi();
const std::vector<std_complex> output_psi();
const std::vector<double> output_n_r();
const complex output_psi(int i);
const std_complex output_psi(int i);
const double output_n_r(int i );
double abs_val();

//FUNKCJA RAPORTUJĄCA O WSZELKICH WARTOŚCIACH NAN
const std::string report_exceptions();
};

0 comments on commit 01d3b49

Please sign in to comment.