Commit 809d1ac3 authored by Alain O' Miniussi's avatar Alain O' Miniussi
Browse files

More assertions


Former-commit-id: b25eab6053a8fbdd0fbf9e015f5023b7927063c5
parent 9b1d9d87
...@@ -50,7 +50,7 @@ namespace D3 { ...@@ -50,7 +50,7 @@ namespace D3 {
PolarGrid *Eta1, PolarGrid *Eta1,
PolarGrid *Eta2, PolarGrid *Eta2,
PolarGrid *Temperature, PolarGrid *Temperature,
PolarGrid *EnergyRad, PolarGrid const& EnergyRad,
PolarGrid *Rhs, PolarGrid *Rhs,
real dt); real dt);
static void ComputeOpticalDepth(PolarGrid *OpaS, static void ComputeOpticalDepth(PolarGrid *OpaS,
...@@ -113,10 +113,9 @@ namespace D3 { ...@@ -113,10 +113,9 @@ namespace D3 {
DiskDesc const& disk = Rho->disk(); DiskDesc const& disk = Rho->disk();
DiskGrid const& grid = disk.grid(); DiskGrid const& grid = disk.grid();
boost::mpi::communicator const& world = grid.comm(); boost::mpi::communicator const& world = grid.comm();
real *temper,* energy, *dens,*energyrad; real *temper,* energy, *dens,*energyrad;
real *eta1,*eta2; real *eta1,*eta2;
static int count=0; static bool firstCall = true;
int nr = Rho->nRad(); int nr = Rho->nRad();
int ns = Rho->nSec(); int ns = Rho->nSec();
...@@ -134,20 +133,24 @@ namespace D3 { ...@@ -134,20 +133,24 @@ namespace D3 {
ComputeEta(Rho,StellarRad,OpaP,Temperature,Viscosity,Eta1,Eta2,psys,dt); ComputeEta(Rho,StellarRad,OpaP,Temperature,Viscosity,Eta1,Eta2,psys,dt);
if (!(Restart || RestartGasOnly)) { if (!(Restart || RestartGasOnly)) {
if (count == 0){ if (firstCall){
firstCall = false;
for (int i = 0; i < nr ; i++) { for (int i = 0; i < nr ; i++) {
for (int h =0; h< ni; h++) { for (int h =0; h< ni; h++) {
for (int j = 0; j < ns; j++) { for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns; int l = j+h*ns+i*ni*ns;
energyrad[l] = (temper[l] - eta1[l])/eta2[l]; energyrad[l] = (temper[l] - eta1[l])/eta2[l];
assert(energyrad[l]>=real(0));
} }
} }
} }
//EnergyRad->values() = (Temperature->values() - Eta1->values())/Eta2->values();
if(grid.first()) { if(grid.first()) {
for (int h =0; h<ni; h++){ for (int h =0; h<ni; h++){
for (int j = 0; j < ns; j++) { for (int j = 0; j < ns; j++) {
int l = j +h*ns; int l = j +h*ns;
energyrad[l]=energyrad[l+ni*ns]; energyrad[l]=energyrad[l+ni*ns];
assert(energyrad[l]>=0);
} }
} }
} }
...@@ -157,6 +160,7 @@ namespace D3 { ...@@ -157,6 +160,7 @@ namespace D3 {
for (int j = 0; j < ns; j++) { for (int j = 0; j < ns; j++) {
int l = j +h*ns+i*ni*ns; int l = j +h*ns+i*ni*ns;
energyrad[l]=energyrad[l-ni*ns]; energyrad[l]=energyrad[l-ni*ns];
assert(energyrad[l]>=0);
} }
} }
} }
...@@ -166,6 +170,7 @@ namespace D3 { ...@@ -166,6 +170,7 @@ namespace D3 {
for (int j = 0; j < ns; j++) { for (int j = 0; j < ns; j++) {
int l = j +h*ns+i*ni*ns; int l = j +h*ns+i*ni*ns;
energyrad[l]= BoundEcode; energyrad[l]= BoundEcode;
assert(energyrad[l]>=0);
} }
} }
} }
...@@ -179,13 +184,15 @@ namespace D3 { ...@@ -179,13 +184,15 @@ namespace D3 {
} else { } else {
energyrad[l]=BoundEcode; energyrad[l]=BoundEcode;
} }
assert(energyrad[l]>=0);
} }
} }
} }
count = 1;
} }
assert(EnergyRad->values().min() >= 0);
} }
ComputeMatrixElements(Rho,OpaR,OpaP,Eta1,Eta2,Temperature,EnergyRad,Rhs,dt); assert(EnergyRad->values().min() >= 0);
ComputeMatrixElements(Rho,OpaR,OpaP,Eta1,Eta2,Temperature,*EnergyRad,Rhs,dt);
SolveMatrix(Rhs,EnergyRad,EnergyRadNew,Temperature,Rho,dt); SolveMatrix(Rhs,EnergyRad,EnergyRadNew,Temperature,Rho,dt);
for (int i = 1; i < nr-1 ; i++) { for (int i = 1; i < nr-1 ; i++) {
...@@ -198,7 +205,6 @@ namespace D3 { ...@@ -198,7 +205,6 @@ namespace D3 {
} }
} }
} }
} }
...@@ -638,7 +644,7 @@ namespace D3 { ...@@ -638,7 +644,7 @@ namespace D3 {
PolarGrid *Eta1, PolarGrid *Eta1,
PolarGrid *Eta2, PolarGrid *Eta2,
PolarGrid *Temperature, PolarGrid *Temperature,
PolarGrid *EnergyRad, PolarGrid const& EnergyRad,
PolarGrid *Rhs, PolarGrid *Rhs,
real dt) real dt)
{ {
...@@ -653,7 +659,7 @@ namespace D3 { ...@@ -653,7 +659,7 @@ namespace D3 {
real *opar = OpaR ->field(); real *opar = OpaR ->field();
real *opap = OpaP ->field(); real *opap = OpaP ->field();
real *temper = Temperature->field(); real *temper = Temperature->field();
real *enrad = EnergyRad->field(); real const *enrad = EnergyRad.field();
real *eta1 = Eta1 ->field(); real *eta1 = Eta1 ->field();
real *eta2 = Eta2 ->field(); real *eta2 = Eta2 ->field();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment