From bc28246bb2a44a3291ffb2440a6395fef07cd151 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Asia=20J=C4=99drzejewska-Szmek?= <asia@in.waw.pl>
Date: Fri, 20 Jan 2017 15:49:32 -0500
Subject: [PATCH] Small changes to make the code more alike genesis and other
 moose objects

Make sure that innerArea is not negative. Move sending out of messages
to the bottom of the process function. Add sending out of
concentration to the buffer and diffusion to other difshells in
DifShell and diffusion to other buffers in DifBuffer in reinit. This
matches genesis implementation.
---
 moose-core/biophysics/DifBuffer.cpp | 29 ++++++++++++++++++-----------
 moose-core/biophysics/DifShell.cpp  | 20 +++++++++++++-------
 2 files changed, 31 insertions(+), 18 deletions(-)

diff --git a/moose-core/biophysics/DifBuffer.cpp b/moose-core/biophysics/DifBuffer.cpp
index da0ebf8b..666bc8d2 100644
--- a/moose-core/biophysics/DifBuffer.cpp
+++ b/moose-core/biophysics/DifBuffer.cpp
@@ -369,9 +369,7 @@ void DifBuffer::vProcess( const Eref & e, ProcPtr p )
    * then compute their incoming fluxes.
    */
 
-  innerDifSourceOut()->send( e, prevFree_, thickness_ );
-  outerDifSourceOut()->send( e, prevFree_, thickness_ );
-  reactionOut()->send(e,kf_,kb_,bFree_,bBound_);
+  
   
   Af_ += kb_ * bBound_;
   Bf_ += kf_ * activation_;
@@ -381,25 +379,31 @@ void DifBuffer::vProcess( const Eref & e, ProcPtr p )
   prevFree_ = bFree_;
   prevBound_ = bBound_;
  
-  Af_ = 0;
-  Bf_= 0;
   /**
    * Send ion concentration to ion buffers. They will send back information on
    * the reaction (forward / backward rates ; free / bound buffer concentration)
    * immediately, which this DifShell will use to find amount of ion captured
    * or released in the current time-step.
    */
+  innerDifSourceOut()->send( e, prevFree_, thickness_ );
+  outerDifSourceOut()->send( e, prevFree_, thickness_ );
+  reactionOut()->send(e,kf_,kb_,bFree_,bBound_);
+  Af_ = 0;
+  Bf_= 0;
 
 }
 
 void DifBuffer::vReinit( const Eref& e, ProcPtr p )
 {
 	
+  Af_ = 0;
+  Bf_= 0;
+  double rOut = diameter_/2.;
   
-  const double rOut = diameter_/2.;
-  
-  const double rIn = rOut - thickness_;
-  
+  double rIn = rOut - thickness_;
+
+  if (rIn<0)
+	 rIn = 0.;
   switch ( shapeMode_ )
     {
       /*
@@ -443,18 +447,21 @@ void DifBuffer::vReinit( const Eref& e, ProcPtr p )
   prevFree_ = bFree_;
   bBound_ = bTot_ - bFree_;
   prevBound_ = bBound_;
+  innerDifSourceOut()->send( e, prevFree_, thickness_ );
+  outerDifSourceOut()->send( e, prevFree_, thickness_ );
+  
 }
 
 void DifBuffer::vFluxFromIn(const Eref& e,double innerC, double innerThickness)
 {
-  double dif = 2 * D_ * innerArea_ / ((thickness_ + innerThickness) * volume_);
+  double dif = 2 * D_ * innerArea_ / (thickness_ + innerThickness)/ volume_;
   Af_ += dif * innerC;
   Bf_ += dif;
 }
 
 void DifBuffer::vFluxFromOut(const Eref& e,double outerC, double outerThickness)
 {
-  double dif = 2 * D_ * outerArea_ / ((thickness_ + outerThickness)  * volume_);
+  double dif = 2 * D_ * outerArea_ / (thickness_ + outerThickness) / volume_;
   Af_ += dif * outerC;
   Bf_ += dif;
 }
diff --git a/moose-core/biophysics/DifShell.cpp b/moose-core/biophysics/DifShell.cpp
index 015f6117..bc893d76 100644
--- a/moose-core/biophysics/DifShell.cpp
+++ b/moose-core/biophysics/DifShell.cpp
@@ -281,9 +281,12 @@ void DifShell::vReinit( const Eref& e, ProcPtr p )
   dCbyDt_ = leak_;
   Cmultiplier_ = 0;
 
-  const double rOut = diameter_/2.;
+  double rOut = diameter_/2.;
   
-  const double rIn = rOut - thickness_;
+  double rIn = rOut - thickness_;
+
+  if (rIn <0)
+	  rIn = 0.;
   
   switch ( shapeMode_ )
     {
@@ -326,6 +329,8 @@ void DifShell::vReinit( const Eref& e, ProcPtr p )
   C_= Ceq_;
   prevC_ = Ceq_;
   concentrationOut()->send( e, C_ );
+  innerDifSourceOut()->send( e, prevC_, thickness_ );
+  outerDifSourceOut()->send( e, prevC_, thickness_ );
 }
 
 void DifShell::vProcess( const Eref & e, ProcPtr p )
@@ -335,9 +340,7 @@ void DifShell::vProcess( const Eref & e, ProcPtr p )
    * then compute their incoming fluxes.
    */
   
-  innerDifSourceOut()->send( e, prevC_, thickness_ );
-  outerDifSourceOut()->send( e, prevC_, thickness_ );
-  concentrationOut()->send( e, C_ );
+ 
   C_ = integrate(C_,p->dt,dCbyDt_,Cmultiplier_);
  
   /**
@@ -351,6 +354,9 @@ void DifShell::vProcess( const Eref & e, ProcPtr p )
   //cout<<"Shell "<< C_<<" ";
   dCbyDt_ = leak_;
   Cmultiplier_ = 0;
+  innerDifSourceOut()->send( e, prevC_, thickness_ );
+  outerDifSourceOut()->send( e, prevC_, thickness_ );
+  concentrationOut()->send( e, C_ );
 
 }
 void DifShell::vBuffer(const Eref& e,
@@ -365,7 +371,7 @@ void DifShell::vBuffer(const Eref& e,
 
 void DifShell::vFluxFromOut(const Eref& e, double outerC, double outerThickness )
 {
-  double diff =2.* ( D_ / volume_ ) * ( outerArea_ / (outerThickness + thickness_) );
+  double diff =2.*  D_  *  outerArea_ / (outerThickness + thickness_) /volume_;
   //influx from outer shell
   /**
    * We could pre-compute ( D / Volume ), but let us leave the optimizations
@@ -380,7 +386,7 @@ void DifShell::vFluxFromIn(const Eref& e, double innerC, double innerThickness )
 {
   //influx from inner shell
   //double dx = ( innerThickness + thickness_ ) / 2.0;
-  double diff = 2.*( D_ / volume_ ) * ( innerArea_ / (innerThickness + thickness_) );
+  double diff = 2.* D_ * innerArea_ / (innerThickness + thickness_) /volume_;
   //cout << "FluxFromIn "<<innerC<<" "<<innerThickness;
   dCbyDt_ +=  diff *  innerC ;
   Cmultiplier_ += diff ;
-- 
GitLab