From f32bef60a5a4f9c0e3747c29347bc500312611c5 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:31 -0500
Subject: [PATCH] Add state in t(i-1)
---
moose-core/biophysics/DifBuffer.cpp | 46 +++++++++++++++++++----------
moose-core/biophysics/DifBuffer.h | 2 ++
moose-core/biophysics/DifShell.cpp | 17 +++++++----
moose-core/biophysics/DifShell.h | 1 +
4 files changed, 46 insertions(+), 20 deletions(-)
diff --git a/moose-core/biophysics/DifBuffer.cpp b/moose-core/biophysics/DifBuffer.cpp
index 0d695f80..5b460fd1 100644
--- a/moose-core/biophysics/DifBuffer.cpp
+++ b/moose-core/biophysics/DifBuffer.cpp
@@ -89,6 +89,8 @@ DifBuffer::DifBuffer() :
Bf_(0),
bFree_(0),
bBound_(0),
+ prevFree_(0),
+ prevBound_(0),
bTot_(0),
kf_(0),
kb_(0),
@@ -137,6 +139,9 @@ void DifBuffer::vSetBFree(const Eref& e,double value)
}
bFree_ = value;
bBound_ = bTot_ - bFree_;
+ prevFree_= bFree_;
+ prevBound_ = bBound_;
+
}
double DifBuffer::vGetBBound(const Eref& e) const
@@ -156,6 +161,8 @@ void DifBuffer::vSetBBound(const Eref& e,double value)
}
bBound_ = value;
bFree_ = bTot_ - bBound_;
+ prevFree_= bFree_;
+ prevBound_ = bBound_;
}
@@ -172,6 +179,7 @@ void DifBuffer::vSetBTot(const Eref& e,double value)
}
bTot_ = value;
bFree_ = bTot_;
+ bBound_ = 0;
}
@@ -340,6 +348,7 @@ void DifBuffer::vBuffer(const Eref& e,
double C )
{
activation_ = C;
+ //cout<<"Buffer"<< C<<" ";
}
@@ -360,14 +369,19 @@ void DifBuffer::vProcess( const Eref & e, ProcPtr p )
* then compute their incoming fluxes.
*/
- innerDifSourceOut()->send( e, bFree_, thickness_ );
- outerDifSourceOut()->send( e, bFree_, thickness_ );
+ innerDifSourceOut()->send( e, prevFree_, thickness_ );
+ outerDifSourceOut()->send( e, prevFree_, thickness_ );
+ reactionOut()->send(e,kf_,kb_,bFree_,bBound_);
+
Af_ += kb_ * bBound_;
Bf_ += kf_ * activation_;
+
bFree_ = integrate(bFree_,p->dt,Af_,Bf_);
bBound_ = bTot_ - bFree_;
-
- reactionOut()->send(e,kf_,kb_,bFree_,bBound_);
+ prevFree_ = bFree_;
+ prevBound_ = bBound_;
+
+
/**
* Send ion concentration to ion buffers. They will send back information on
* the reaction (forward / backward rates ; free / bound buffer concentration)
@@ -380,10 +394,11 @@ void DifBuffer::vProcess( const Eref & e, ProcPtr p )
void DifBuffer::vReinit( const Eref& e, ProcPtr p )
{
- const double dOut = diameter_;
- const double dIn = diameter_ - thickness_;
- bFree_ = bTot_;
- bBound_ = 0;
+
+ const double rOut = diameter_/2.;
+
+ const double rIn = rOut - thickness_;
+
switch ( shapeMode_ )
{
/*
@@ -391,15 +406,15 @@ void DifBuffer::vReinit( const Eref& e, ProcPtr p )
*/
case 0:
if ( length_ == 0.0 ) { // Spherical shell
- volume_ = ( M_PI / 6.0 ) * ( dOut * dOut * dOut - dIn * dIn * dIn );
- outerArea_ = M_PI * dOut * dOut;
- innerArea_ = M_PI * dIn * dIn;
+ volume_ = 4./3.* M_PI * ( rOut * rOut * rOut - rIn * rIn * rIn );
+ outerArea_ = 4*M_PI * rOut * rOut;
+ innerArea_ = 4*M_PI * rIn * rIn;
} else { // Cylindrical shell
- volume_ = ( M_PI * length_ / 4.0 ) * ( dOut * dOut - dIn * dIn );
- outerArea_ = M_PI * dOut * length_;
- innerArea_ = M_PI * dIn * length_;
+ volume_ = ( M_PI * length_ ) * ( rOut * rOut - rIn * rIn );
+ outerArea_ = 2*M_PI * rOut * length_;
+ innerArea_ = 2*M_PI * rIn * length_;
}
-
+
break;
/*
@@ -422,6 +437,7 @@ void DifBuffer::vReinit( const Eref& e, ProcPtr p )
default:
assert( 0 );
}
+
}
void DifBuffer::vFluxFromIn(const Eref& e,double innerC, double innerThickness)
diff --git a/moose-core/biophysics/DifBuffer.h b/moose-core/biophysics/DifBuffer.h
index fe3f7772..60f57445 100644
--- a/moose-core/biophysics/DifBuffer.h
+++ b/moose-core/biophysics/DifBuffer.h
@@ -74,6 +74,8 @@ class DifBuffer: public DifBufferBase{
double Bf_;
double bFree_; //free buffer concentration
double bBound_; //bound buffer concentration
+ double prevFree_;
+ double prevBound_;
//double prevFree_; //bFree at previous time
//double prevBound_; //bBound at previous time
double bTot_; //bFree+bBound
diff --git a/moose-core/biophysics/DifShell.cpp b/moose-core/biophysics/DifShell.cpp
index c0e3040b..d4e37d08 100644
--- a/moose-core/biophysics/DifShell.cpp
+++ b/moose-core/biophysics/DifShell.cpp
@@ -58,6 +58,7 @@ DifShell::DifShell() :
Cmultiplier_(0.0),
C_( 0.0 ),
Ceq_( 0.0 ),
+ prevC_(0.0),
D_( 0.0 ),
valence_( 0.0 ),
leak_( 0.0 ),
@@ -84,6 +85,7 @@ void DifShell::vSetC( const Eref& e, double C)
}
C_ = C;
+ prevC_ = C_;
}
double DifShell::vGetC(const Eref& e) const
{
@@ -270,6 +272,7 @@ double DifShell::integrate( double state, double dt, double A, double B )
double x = exp( -B * dt );
return state * x + ( A / B ) * ( 1 - x );
}
+
return state + A * dt ;
}
@@ -277,7 +280,7 @@ void DifShell::vReinit( const Eref& e, ProcPtr p )
{
dCbyDt_ = leak_;
Cmultiplier_ = 0;
- C_ = Ceq_;
+
const double rOut = diameter_/2.;
const double rIn = rOut - thickness_;
@@ -329,19 +332,23 @@ void DifShell::vProcess( const Eref & e, ProcPtr p )
* then compute their incoming fluxes.
*/
- innerDifSourceOut()->send( e, C_, thickness_ );
- outerDifSourceOut()->send( e, C_, thickness_ );
+ innerDifSourceOut()->send( e, prevC_, thickness_ );
+ outerDifSourceOut()->send( e, prevC_, thickness_ );
+ concentrationOut()->send( e, C_ );
C_ = integrate(C_,p->dt,dCbyDt_,Cmultiplier_);
-
+
/**
* 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.
*/
- concentrationOut()->send( e, C_ );
+ prevC_ = C_;
+
+ //cout<<"Shell "<< C_<<" ";
dCbyDt_ = leak_;
Cmultiplier_ = 0;
+
}
void DifShell::vBuffer(const Eref& e,
double kf,
diff --git a/moose-core/biophysics/DifShell.h b/moose-core/biophysics/DifShell.h
index 6bef0803..42cf643d 100644
--- a/moose-core/biophysics/DifShell.h
+++ b/moose-core/biophysics/DifShell.h
@@ -84,6 +84,7 @@ class DifShell: public DifShellBase{
double Cmultiplier_;
double C_;
double Ceq_;
+ double prevC_;
double D_;
double valence_;
double leak_;
--
GitLab