[cig-commits] [commit] master: Moved many variables from Block class to Sim (48b16b5)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 2 13:56:07 PST 2014
Repository : https://github.com/geodynamics/vq
On branch : master
Link : https://github.com/geodynamics/vq/compare/06902a05cf4d70bc94c85d1c195c073d4e74cdad...c935dfd33f870a081b6c01cce97a060ad1cbdbda
>---------------------------------------------------------------
commit 48b16b5be41f4ed85754090ff23590e08c1cfeea
Author: Eric Heien <emheien at ucdavis.edu>
Date: Tue Nov 25 12:37:02 2014 -0800
Moved many variables from Block class to Sim
>---------------------------------------------------------------
48b16b5be41f4ed85754090ff23590e08c1cfeea
quakelib/src/QuakeLib.h | 2 +
src/core/Block.cpp | 57 +--------------
src/core/Block.h | 131 -----------------------------------
src/core/InitBlocks.cpp | 7 --
src/core/SimData.cpp | 55 +++++++++++++--
src/core/SimData.h | 54 ++++++---------
src/core/SimDataBlocks.cpp | 2 -
src/core/Simulation.cpp | 34 +++++----
src/core/Simulation.h | 110 ++++++++++++++++++++++++++++-
src/io/CheckpointFileOutput.cpp | 2 +-
src/io/HDF5Data.cpp | 20 +++---
src/io/ReadModelFile.cpp | 7 +-
src/simulation/GreensInit.cpp | 2 +
src/simulation/ProgressMonitor.cpp | 2 +-
src/simulation/RunEvent.cpp | 65 +++++++++--------
src/simulation/SanityCheck.cpp | 2 +-
src/simulation/UpdateBlockStress.cpp | 32 ++++++---
17 files changed, 267 insertions(+), 317 deletions(-)
diff --git a/quakelib/src/QuakeLib.h b/quakelib/src/QuakeLib.h
index 94fba20..e74b072 100644
--- a/quakelib/src/QuakeLib.h
+++ b/quakelib/src/QuakeLib.h
@@ -129,6 +129,7 @@ namespace quakelib {
//! Set the slip rate in m/s for this block.
void set_slip_rate(const double &new_slip_rate) throw(std::invalid_argument) {
if (isnan(new_slip_rate)) throw std::invalid_argument("quakelib::Element::set_slip_rate");
+
_slip_rate = new_slip_rate;
};
//! Get the slip rate in cm/year for this block.
@@ -139,6 +140,7 @@ namespace quakelib {
//! Set the rake angle of this block in radians.
void set_rake(const double &new_rake) throw(std::invalid_argument) {
if (isnan(new_rake)) throw std::invalid_argument("quakelib::Element::set_rake");
+
_rake = new_rake;
};
//! Get the rake angle of this block in radians.
diff --git a/src/core/Block.cpp b/src/core/Block.cpp
index f66994a..848317a 100644
--- a/src/core/Block.cpp
+++ b/src/core/Block.cpp
@@ -126,60 +126,5 @@ void Block::clear(void) {
id = UNDEFINED_ELEMENT_ID;
fid = UNDEFINED_FAULT_ID;
sid = UNDEFINED_SECTION_ID;
- self_shear = self_normal = dynamic_val = stress_drop = std::numeric_limits<float>::quiet_NaN();
- init_shear_stress = init_normal_stress = rhogd = friction_val = std::numeric_limits<float>::quiet_NaN();
-}
-
-void Block::calcFriction(void) {
- friction_val = (fabs(getStressDrop())/getRhogd());
-}
-
-bool Block::cffFailure(void) const {
- return (state.cff >= 0);
-}
-
-bool Block::dynamicFailure(const FaultID &event_fault) const {
- return (fid == event_fault &&
- state.cff > state.cff0 &&
- //state.cff > getStressDrop() &&
- fabs((state.cff0-state.cff)/state.cff0) > dynamic_val);
-}
-
-void Block::calcCFF(void) {
- state.cff = fabs(state.shear_stress[0]) - fabs(friction()*state.normal_stress[0]);
-}
-
-void Block::setStatePtrs(double *shear_stress, double *normal_stress, double *update_field) {
- state.shear_stress = shear_stress;
- state.normal_stress = normal_stress;
- state.updateField = update_field;
-}
-
-StateCheckpointData State::readCheckpointData(void) const {
- StateCheckpointData ckpt_data;
-
- ckpt_data.slipDeficit = slipDeficit;
- ckpt_data.cff = cff;
- ckpt_data.shear_stress = *shear_stress;
- ckpt_data.normal_stress = *normal_stress;
- ckpt_data.updateField = *updateField;
-
- return ckpt_data;
-}
-
-void State::storeCheckpointData(const StateCheckpointData &ckpt_data) {
- slipDeficit = ckpt_data.slipDeficit;
- cff = ckpt_data.cff;
- *shear_stress = ckpt_data.shear_stress;
- *normal_stress = ckpt_data.normal_stress;
- *updateField = ckpt_data.updateField;
-}
-
-std::ostream &operator<<(std::ostream &os, const State &s) {
- os << std::setw(10) << s.cff
- << std::setw(12) << s.updateField[0]
- << std::setw(10) << s.shear_stress[0]
- << std::setw(10) << s.normal_stress[0];
-
- return os;
+ dynamic_val = init_shear_stress = init_normal_stress = std::numeric_limits<float>::quiet_NaN();
}
diff --git a/src/core/Block.h b/src/core/Block.h
index 4ec370a..9df2b52 100644
--- a/src/core/Block.h
+++ b/src/core/Block.h
@@ -72,36 +72,6 @@ typedef struct BlockSweepVals BlockSweepVals;
std::ostream &operator<<(std::ostream &os, const BlockVal &bv);
/*!
- Class to keep dynamic block attributes. Not all of them are needed to save between simulations,
- only slipDeficit, but they are saved anyway to provide debugging and other quickly accessible info.
- */
-class State {
- public:
- State(void) : slipDeficit(0), cff(0), cff0(0), stressS0(0), stressN0(0), shear_stress(NULL), normal_stress(NULL), updateField(NULL) {};
- //! slip - Vt
- double slipDeficit;
- //! coulomb failure function
- double cff;
- //! cff before an event
- double cff0;
- //! shear stress before an event
- double stressS0;
- //! normal stress before an event
- double stressN0;
- //! ptr to shear stress in simulation
- double *shear_stress;
- //! ptr to normal stress in simulation
- double *normal_stress;
- //! ptr to update_field in simulation
- double *updateField;
-
- StateCheckpointData readCheckpointData(void) const;
- void storeCheckpointData(const StateCheckpointData &ckpt_data);
-
- friend std::ostream &operator<<(std::ostream &os, const State &b);
-};
-
-/*!
Class of block attributes. These attributes are static during the simulation except rand.
*/
class Block : public quakelib::SimElement {
@@ -111,45 +81,20 @@ class Block : public quakelib::SimElement {
SectionID sid; // section id
// constants during simulation
- double self_shear; // Greens function shear by block on itself
- double self_normal; // Greens function normal by block on itself
double dynamic_val; // dynamic failure value for this block
- double stress_drop; // predetermined stress drop (in pascals)
-
double init_shear_stress; // initial shear stress
double init_normal_stress; // initial normal stress
- double rhogd; // normal stress
-
- double friction_val; // friction value for this block
-
bool failed;
public:
void get_rake_and_normal_stress_due_to_block(double stresses[2], const double &sample_dist, const Block &source_block) const;
- State state;
-
//! Resets the values of this block.
void clear(void);
- //! Calculates the static friction of this block.
- void calcFriction(void);
- //! Get the static friction of this block.
- double friction(void) const {
- return friction_val;
- };
-
- //! Whether the block experienced static friction failure.
- //! This occurs if the Coulomb failure function goes over 0.
- bool cffFailure(void) const;
- //! Whether the block experienced dynamic friction failure.
- //! This occurs if the block is on the same fault as the trigger
- //! and it undergoes a significant CFF shift during this event
- //! (where significant is defined in terms of dynamic_val).
- bool dynamicFailure(const FaultID &event_fault) const;
bool getFailed(void) const {
return failed;
};
@@ -157,45 +102,6 @@ class Block : public quakelib::SimElement {
failed = in_failed;
};
- //! Returns the precalculated CFF of this block.
- double getCFF(void) const {
- return state.cff;
- };
- //! Returns the slip deficit of this block.
- double getSlipDeficit(void) const {
- return state.slipDeficit;
- };
- //! Returns the slip deficit of this block.
- double getShearStress(void) const {
- return state.shear_stress[0];
- };
- //! Returns the slip deficit of this block.
- double getNormalStress(void) const {
- return state.normal_stress[0];
- };
- //! Calculates and stores the CFF of this block.
- void calcCFF(void);
- //! Record the current stresses and CFF.
- void saveStresses(void) {
- state.stressS0 = state.shear_stress[0];
- state.stressN0 = state.normal_stress[0];
- state.cff0 = state.cff;
- };
- //! Get the recorded shear stress.
- double getStressS0(void) const {
- return state.stressS0;
- };
- //! Get the recorded normal stress.
- double getStressN0(void) const {
- return state.stressN0;
- };
- //! Get the recorded CFF.
- double getCFF0(void) const {
- return state.cff0;
- };
- //! Set pointers into arrays for the stress values of this block
- void setStatePtrs(double *shear_stress, double *normal_stress, double *update_field);
-
// Functions for manipulation of block parameters
//! Set the block ID of this block.
void setBlockID(const BlockID &new_id) {
@@ -224,21 +130,6 @@ class Block : public quakelib::SimElement {
return sid;
};
- //! Calculate the expected recurrence of this block in years.
- // TODO: fix this, it is currently incorrect for multiple block models
- double getRecurrence(void) const {
- return stress_drop/(_slip_rate*(self_shear-self_normal*friction_val));
- };
- //! Set the Greens function self shear of this block.
- void setSelfStresses(const double &new_self_shear, const double &new_self_normal) {
- self_shear = new_self_shear;
- self_normal = new_self_normal;
- };
- //! Get the Greens function self shear of this block.
- double getSelfStresses(void) const {
- return self_shear - friction()*self_normal;
- };
-
//! Get the dynamic triggering value for this block.
double getDynamicVal(void) const {
return dynamic_val;
@@ -264,28 +155,6 @@ class Block : public quakelib::SimElement {
double getInitNormalStress(void) const {
return init_normal_stress;
};
-
- //! Get the stress drop for this block in bars.
- double getStressDrop(void) const {
- return stress_drop;
- };
- //! Set the stress drop for this block in bars.
- void setStressDrop(const double &new_stress_drop) {
- stress_drop = new_stress_drop;
- calcFriction();
- };
- double getRhogd(void) const {
- return rhogd;
- };
- void setRhogd(const double &new_rhogd) {
- rhogd = new_rhogd;
- calcFriction();
- };
-
- //! Get the string header for state files.
- static std::string getStateHeader(void) {
- return " cff updateField shear_stress normal_stress";
- };
};
typedef std::vector<Block> BlockList;
diff --git a/src/core/InitBlocks.cpp b/src/core/InitBlocks.cpp
index 49d8b4d..eab82af 100644
--- a/src/core/InitBlocks.cpp
+++ b/src/core/InitBlocks.cpp
@@ -41,7 +41,6 @@ void VCInitBlocks::dryRun(SimFramework *_sim) {
*/
void VCInitBlocks::init(SimFramework *_sim) {
Simulation *sim = static_cast<Simulation *>(_sim);
- int i;
BlockList::iterator bit;
assertThrow(sim->numGlobalBlocks() > 0, "Simulation must include at least 1 block.");
@@ -56,12 +55,6 @@ void VCInitBlocks::init(SimFramework *_sim) {
// transposed array for faster sweep calculations
sim->useTransposedMatrix());
- // Set the block pointers to the storage arrays
- for (i=0,bit=sim->begin(); bit!=sim->end(); ++bit,++i)
- bit->setStatePtrs(sim->getShearStressPtr(i),
- sim->getNormalStressPtr(i),
- sim->getUpdateFieldPtr(i));
-
// Set the starting year of the simulation
sim->setYear(sim->getSimStart());
diff --git a/src/core/SimData.cpp b/src/core/SimData.cpp
index 2a925e5..d9cb7ef 100644
--- a/src/core/SimData.cpp
+++ b/src/core/SimData.cpp
@@ -30,8 +30,6 @@ void VCSimData::setupArrays(const unsigned int &global_sys_size,
const unsigned int &local_sys_size,
const bool &compressed,
const bool &transposed) {
- int i;
-
deallocateArrays();
global_size = global_sys_size;
@@ -62,10 +60,33 @@ void VCSimData::setupArrays(const unsigned int &global_sys_size,
assertThrow(normal_stress, "Not enough memory to allocate normal stress array.");
update_field = (double *)malloc(sizeof(double)*global_size);
assertThrow(update_field, "Not enough memory to allocate update field array.");
-
- for (i=0; i<global_size; ++i) {
- shear_stress[i] = normal_stress[i] = 0;
- update_field[i] = 0;
+ slip_deficit = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(slip_deficit, "Not enough memory to allocate slip deficit array.");
+ rhogd = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(rhogd, "Not enough memory to allocate rhogd array.");
+ stress_drop = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(stress_drop, "Not enough memory to allocate stress drop array.");
+ cff = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(cff, "Not enough memory to allocate cff array.");
+ friction = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(friction, "Not enough memory to allocate friction array.");
+ cff0 = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(cff0, "Not enough memory to allocate cff0 array.");
+ self_shear = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(self_shear, "Not enough memory to allocate self_shear array.");
+ self_normal = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(self_normal, "Not enough memory to allocate self_normal array.");
+ shear_stress0 = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(shear_stress0, "Not enough memory to allocate shear_stress0 array.");
+ normal_stress0 = (double *)malloc(sizeof(double)*global_size);
+ assertThrow(normal_stress0, "Not enough memory to allocate normal_stress0 array.");
+
+ for (BlockID i=0; i<global_size; ++i) {
+ shear_stress[i] = normal_stress[i] = std::numeric_limits<float>::quiet_NaN();
+ update_field[i] = slip_deficit[i] = std::numeric_limits<float>::quiet_NaN();
+ rhogd[i] = stress_drop[i] = cff[i] = std::numeric_limits<float>::quiet_NaN();
+ friction[i] = cff0[i] = self_shear[i] = self_normal[i] = std::numeric_limits<float>::quiet_NaN();
+ shear_stress0[i] = normal_stress0[i] = std::numeric_limits<float>::quiet_NaN();
}
}
@@ -97,6 +118,28 @@ void VCSimData::deallocateArrays(void) {
if (update_field) free(update_field);
+ if (slip_deficit) free(slip_deficit);
+
+ if (rhogd) free(rhogd);
+
+ if (stress_drop) free(stress_drop);
+
+ if (cff) free(cff);
+
+ if (friction) free(friction);
+
+ if (cff0) free(cff0);
+
+ if (self_shear) free(self_shear);
+
+ if (self_normal) free(self_normal);
+
+ if (shear_stress0) free(shear_stress0);
+
+ if (normal_stress0) free(normal_stress0);
+
green_shear = green_normal = NULL;
shear_stress = normal_stress = update_field = NULL;
+ slip_deficit = rhogd = stress_drop = cff = friction = NULL;
+ cff0 = self_shear = self_normal = shear_stress0 = normal_stress0 = NULL;
}
diff --git a/src/core/SimData.h b/src/core/SimData.h
index 520da67..899b1a1 100644
--- a/src/core/SimData.h
+++ b/src/core/SimData.h
@@ -40,13 +40,26 @@ class VCSimData : public VCSimDataBlocks, public VCSimDataEvents {
unsigned int global_size, local_size;
quakelib::DenseMatrix<GREEN_VAL> *green_shear, *green_normal;
+ protected:
double *shear_stress;
double *normal_stress;
double *update_field;
+ double *slip_deficit;
+ double *rhogd;
+ double *stress_drop;
+ double *cff;
+ double *friction;
+ double *cff0;
+ double *self_shear;
+ double *self_normal;
+ double *shear_stress0;
+ double *normal_stress0;
public:
VCSimData(void) : global_size(0), local_size(0), green_shear(NULL), green_normal(NULL),
- shear_stress(NULL), normal_stress(NULL), update_field(NULL) {};
+ shear_stress(NULL), normal_stress(NULL), update_field(NULL), slip_deficit(NULL),
+ rhogd(NULL), stress_drop(NULL), cff(NULL), friction(NULL), cff0(NULL),
+ self_shear(NULL), self_normal(NULL), shear_stress0(NULL), normal_stress0(NULL) {};
void setupArrays(const unsigned int &global_sys_size,
const unsigned int &local_sys_size,
@@ -67,41 +80,14 @@ class VCSimData : public VCSimDataBlocks, public VCSimDataEvents {
quakelib::DenseMatrix<GREEN_VAL> *greenNormal(void) const {
return green_normal;
};
- double *getUpdateFieldPtr(const unsigned int &elem=0) const {
- return &(update_field[elem]);
+ double *getUpdateFieldPtr(void) const {
+ return update_field;
};
- double *getShearStressPtr(const unsigned int &elem=0) const {
- return &(shear_stress[elem]);
+ double *getShearStressPtr(void) const {
+ return shear_stress;
};
- double *getNormalStressPtr(const unsigned int &elem=0) const {
- return &(normal_stress[elem]);
- };
-
- double getUpdateField(const BlockID &b) const {
- return update_field[b];
- };
- void setUpdateField(const BlockID &b, const double &new_val) {
- update_field[b] = new_val;
- };
-
- double getNormalStress(const BlockID &b) const {
- return normal_stress[b];
- };
- void setNormalStress(const BlockID &b, const double &new_val) {
- normal_stress[b] = new_val;
- };
- void addToNormalStress(const BlockID &b, const double &add_val) {
- normal_stress[b] += add_val;
- };
-
- double getShearStress(const BlockID &b) const {
- return shear_stress[b];
- };
- void setShearStress(const BlockID &b, const double &new_val) {
- shear_stress[b] = new_val;
- };
- void addToShearStress(const BlockID &b, const double &add_val) {
- shear_stress[b] += add_val;
+ double *getNormalStressPtr(void) const {
+ return normal_stress;
};
};
diff --git a/src/core/SimDataBlocks.cpp b/src/core/SimDataBlocks.cpp
index 608e88a..d0db699 100644
--- a/src/core/SimDataBlocks.cpp
+++ b/src/core/SimDataBlocks.cpp
@@ -29,12 +29,10 @@ BlockID VCSimDataBlocks::addBlock(const Block &new_block) {
// Perform block sanity checks
// Blocks may not have a negative slip rate (instead, use rake)
- //assertThrow(new_block.dip()>=0&&new_block.dip()<=M_PI/2.0, "Block dip must be between 0 and 90 degrees.");
assertThrow(new_block.getFaultID()>=0, "Block fault ID must be non-negative.");
assertThrow(new_block.min_depth()>new_block.max_depth(), "Block top must be higher than block bottom.");
assertThrow(new_block.slip_rate()>=0, "Blocks may not have a negative slip rate (use rake to specify direction instead).");
assertThrow(new_block.area()>0, "Blocks must have a positive non-zero area.");
- assertThrow(new_block.getRhogd()>0, "Blocks must have a positive rhogd.");
new_block_id = blocks.size();
blocks.push_back(new_block);
diff --git a/src/core/Simulation.cpp b/src/core/Simulation.cpp
index 4922615..403eb65 100644
--- a/src/core/Simulation.cpp
+++ b/src/core/Simulation.cpp
@@ -124,28 +124,28 @@ void Simulation::getInitialFinalStresses(const quakelib::ElementIDSet &block_set
// Add the before/after stresses if it is on this node
// Non-local blocks will have incorrect stress data
if (isLocalToNode(*it)) {
- shear_init += getBlock(*it).getStressS0();
- shear_final += getBlock(*it).getShearStress();
- normal_init += getBlock(*it).getStressN0();
- normal_final += getBlock(*it).getNormalStress();
+ shear_init += shear_stress0[*it];
+ shear_final += getShearStress(*it);
+ normal_init += normal_stress0[*it];
+ normal_final += getNormalStress(*it);
}
}
}
void Simulation::sumStresses(const quakelib::ElementIDSet &block_set,
- double &shear_stress,
- double &shear_stress0,
- double &normal_stress,
- double &normal_stress0) const {
+ double &shear_stress_sum,
+ double &shear_stress0_sum,
+ double &normal_stress_sum,
+ double &normal_stress0_sum) const {
quakelib::ElementIDSet::const_iterator it;
- shear_stress = shear_stress0 = normal_stress = normal_stress0 = 0;
+ shear_stress_sum = shear_stress0_sum = normal_stress_sum = normal_stress0_sum = 0;
for (it=block_set.begin(); it!=block_set.end(); ++it) {
- shear_stress += getShearStress(*it);
- shear_stress0 += getBlock(*it).getStressS0();
- normal_stress += getNormalStress(*it);
- normal_stress0 += getBlock(*it).getStressN0();
+ shear_stress_sum += getShearStress(*it);
+ shear_stress0_sum += shear_stress0[*it];
+ normal_stress_sum += getNormalStress(*it);
+ normal_stress0_sum += normal_stress0[*it];
}
}
@@ -192,10 +192,16 @@ void Simulation::computeCFFs(void) {
int i;
for (i=0; i<numLocalBlocks(); ++i) {
- getBlock(getGlobalBID(i)).calcCFF();
+ BlockID gid = getGlobalBID(i);
+ calcCFF(gid);
}
}
+//! Calculates and stores the CFF of this block.
+void Simulation::calcCFF(const BlockID gid) {
+ cff[gid] = fabs(shear_stress[gid]) - fabs(friction[gid]*normal_stress[gid]);
+}
+
/*!
Performs a matrix-vector multiply (C += A * B) where C is Nx1, A is MxN and B is Nx1.
This directly accesses the matrix A and assumes it is stored in transpose format.
diff --git a/src/core/Simulation.h b/src/core/Simulation.h
index a0efdaa..03828fd 100644
--- a/src/core/Simulation.h
+++ b/src/core/Simulation.h
@@ -77,7 +77,7 @@ class Simulation : public SimFramework, public VCParams, public VCSimData, publi
int numFaults(void) const;
void getInitialFinalStresses(const quakelib::ElementIDSet &block_set, double &shear_init, double &shear_final, double &normal_init, double &normal_final) const;
- void sumStresses(const quakelib::ElementIDSet &block_set, double &shear_stress, double &shear_stress0, double &normal_stress, double &normal_stress0) const;
+ void sumStresses(const quakelib::ElementIDSet &block_set, double &shear_stress_sum, double &shear_stress0_sum, double &normal_stress_sum, double &normal_stress0_sum) const;
double sumGreenShear(const BlockID &r) const {
double sum = 0;
@@ -93,7 +93,7 @@ class Simulation : public SimFramework, public VCParams, public VCSimData, publi
greenShear()->setVal(getLocalInd(r), c, new_green_shear);
greenNormal()->setVal(getLocalInd(r), c, new_green_normal);
- if (r == c) getBlock(r).setSelfStresses(new_green_shear, new_green_normal);
+ if (r == c) setSelfStresses(r, new_green_shear, new_green_normal);
};
double getGreenNormal(const BlockID &r, const BlockID &c) const {
return greenNormal()->val(getLocalInd(r), c);
@@ -116,8 +116,8 @@ class Simulation : public SimFramework, public VCParams, public VCSimData, publi
};
void determineBlockNeighbors(void);
-
void computeCFFs(void);
+ void calcCFF(const BlockID gid);
void matrixVectorMultiplyAccum(double *c, const quakelib::DenseMatrix<GREEN_VAL> *a, const double *b, const bool dense);
void multiplySumRow(double *c, const double *b, const GREEN_VAL *a, const int n, const bool dense);
void multiplyRow(double *c, const double *b, const GREEN_VAL *a, const int n);
@@ -132,6 +132,110 @@ class Simulation : public SimFramework, public VCParams, public VCSimData, publi
return (block_node_map.at(block_id) == node_rank);
};
+ double getSlipDeficit(const BlockID gid) {
+ return slip_deficit[gid];
+ };
+ void setSlipDeficit(const BlockID gid, const double new_slip_deficit) {
+ slip_deficit[gid] = new_slip_deficit;
+ };
+
+ double getCFF(const BlockID gid) {
+ return cff[gid];
+ };
+ void setCFF(const BlockID gid, const double new_cff) {
+ cff[gid] = new_cff;
+ };
+ double getCFF0(const BlockID gid) {
+ return cff0[gid];
+ };
+
+ //! Calculates the static friction of this block.
+ void calcFriction(const BlockID gid) {
+ friction[gid] = (fabs(stress_drop[gid])/rhogd[gid]);
+ }
+ //! Get the static friction of this block.
+ double getFriction(const BlockID gid) {
+ return friction[gid];
+ }
+
+ //! Whether the block experienced static friction failure.
+ //! This occurs if the Coulomb failure function goes over 0.
+ bool cffFailure(const BlockID gid) const {
+ return (cff[gid] >= 0);
+ }
+
+ //! Whether the block experienced dynamic friction failure.
+ //! This occurs if the block is on the same fault as the trigger
+ //! and it undergoes a significant CFF shift during this event
+ //! (where significant is defined in terms of dynamic_val).
+ bool dynamicFailure(const BlockID gid, const FaultID event_fault) const {
+ return (getBlock(gid).getFaultID() == event_fault &&
+ cff[gid] > cff0[gid] &&
+ //state.cff > getStressDrop() &&
+ fabs((cff0[gid]-cff[gid])/cff0[gid]) > getBlock(gid).getDynamicVal());
+ }
+
+ //! Calculate the expected recurrence of this block in years.
+ // TODO: fix this, it is currently incorrect for multiple block models
+ double getRecurrence(const BlockID gid) const {
+ return stress_drop[gid]/(getBlock(gid).slip_rate()*(self_shear[gid]-self_normal[gid]*friction[gid]));
+ };
+
+ //! Set the Greens function self shear of this block.
+ void setSelfStresses(const BlockID gid, const double new_self_shear, const double new_self_normal) {
+ self_shear[gid] = new_self_shear;
+ self_normal[gid] = new_self_normal;
+ };
+ //! Get the Greens function self shear of this block.
+ double getSelfStresses(const BlockID gid) const {
+ return self_shear[gid] - friction[gid]*self_normal[gid];
+ };
+
+ //! Get the stress drop for this block in bars.
+ double getStressDrop(const BlockID gid) const {
+ return stress_drop[gid];
+ };
+ //! Set the stress drop for this block in bars.
+ void setStressDrop(const BlockID gid, const double new_stress_drop) {
+ stress_drop[gid] = new_stress_drop;
+ calcFriction(gid);
+ };
+
+ double getRhogd(const BlockID gid) {
+ return rhogd[gid];
+ };
+ void setRhogd(const BlockID gid, const double new_rhogd) {
+ rhogd[gid] = new_rhogd;
+ calcFriction(gid);
+ };
+
+ double getUpdateField(const BlockID &b) const {
+ return update_field[b];
+ };
+ void setUpdateField(const BlockID &b, const double new_val) {
+ update_field[b] = new_val;
+ };
+
+ double getNormalStress(const BlockID &b) const {
+ return normal_stress[b];
+ };
+ void setNormalStress(const BlockID &b, const double &new_val) {
+ normal_stress[b] = new_val;
+ };
+
+ double getShearStress(const BlockID &b) const {
+ return shear_stress[b];
+ };
+ void setShearStress(const BlockID &b, const double &new_val) {
+ shear_stress[b] = new_val;
+ };
+
+ void saveStresses(const BlockID bid) {
+ shear_stress0[bid] = shear_stress[bid];
+ normal_stress0[bid] = normal_stress[bid];
+ cff0[bid] = cff[bid];
+ }
+
void partitionBlocks(void);
private:
diff --git a/src/io/CheckpointFileOutput.cpp b/src/io/CheckpointFileOutput.cpp
index 2ab56be..9e395a4 100644
--- a/src/io/CheckpointFileOutput.cpp
+++ b/src/io/CheckpointFileOutput.cpp
@@ -35,7 +35,7 @@ void CheckpointFileOutput::writeCheckpoint(const std::string &ckpt_file_name, co
// Get current state for local blocks
for (i=0; i<sim->numLocalBlocks(); ++i) {
bid = sim->getGlobalBID(i);
- checkpoint_set[bid] = sim->getBlock(bid).state.readCheckpointData();
+ //checkpoint_set[bid] = sim->getBlock(bid).state.readCheckpointData();
}
#ifdef HDF5_FOUND
diff --git a/src/io/HDF5Data.cpp b/src/io/HDF5Data.cpp
index 054d461..4d3e1d3 100644
--- a/src/io/HDF5Data.cpp
+++ b/src/io/HDF5Data.cpp
@@ -59,7 +59,7 @@ HDF5CheckpointWriter::HDF5CheckpointWriter(const std::string &ckpt_file_name,
const double &cur_year,
const unsigned int &cur_event,
const CheckpointSet &checkpoints) : HDF5Checkpoint() {
- hsize_t dims[2] = {nblocks, CHECKPOINT_NUM_ENTRIES_HDF5};
+ /*hsize_t dims[2] = {nblocks, CHECKPOINT_NUM_ENTRIES_HDF5};
hsize_t single_val[1] = {1};
herr_t status;
CheckpointSet::const_iterator it;
@@ -74,11 +74,11 @@ HDF5CheckpointWriter::HDF5CheckpointWriter(const std::string &ckpt_file_name,
if (plist_id < 0) exit(-1);
-#ifdef MPI_C_FOUND
-#ifdef H5_HAVE_PARALLEL
+ #ifdef MPI_C_FOUND
+ #ifdef H5_HAVE_PARALLEL
H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL);
-#endif
-#endif
+ #endif
+ #endif
// Create the data file, overwriting any old files
data_file = H5Fcreate(ckpt_file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
@@ -138,11 +138,11 @@ HDF5CheckpointWriter::HDF5CheckpointWriter(const std::string &ckpt_file_name,
// Write all block state data in parallel
xfer_plist_id = H5Pcreate(H5P_DATASET_XFER);
-#ifdef MPI_C_FOUND
-#ifdef H5_HAVE_PARALLEL
+ #ifdef MPI_C_FOUND
+ #ifdef H5_HAVE_PARALLEL
H5Pset_dxpl_mpio(xfer_plist_id, H5FD_MPIO_COLLECTIVE);
-#endif
-#endif
+ #endif
+ #endif
status = H5Dwrite(state_dataset, H5T_NATIVE_DOUBLE, mem_select, file_select, xfer_plist_id, mem_state);
res = H5Sclose(file_select);
@@ -159,7 +159,7 @@ HDF5CheckpointWriter::HDF5CheckpointWriter(const std::string &ckpt_file_name,
res = H5Fclose(data_file);
- delete mem_state;
+ delete mem_state;*/
}
/*!
diff --git a/src/io/ReadModelFile.cpp b/src/io/ReadModelFile.cpp
index 6e1d60b..42e8936 100644
--- a/src/io/ReadModelFile.cpp
+++ b/src/io/ReadModelFile.cpp
@@ -89,15 +89,10 @@ void ReadModelFile::init(SimFramework *_sim) {
new_block.setSectionID(eit->section_id()); // TODO: add sections?
new_block.setFailed(false);
- double rho = 5.515e3; // density of rock in kg m^-3
- double g = 9.81; // force of gravity in m s^-2
- double depth = -new_block.center()[2]; // depth of block center in m
-
- new_block.setRhogd(rho*g*depth); // kg m^-3 * m s^-2 * m = kg m^-1 * s^-2 = Pa
new_block.setInitShearStress(0);
new_block.setInitNormalStress(0);
new_block.setDynamicVal(sim->getDynamic());
- sim->addBlock(new_block);
+ BlockID new_id = sim->addBlock(new_block);
}
}
diff --git a/src/simulation/GreensInit.cpp b/src/simulation/GreensInit.cpp
index 794d894..989e762 100644
--- a/src/simulation/GreensInit.cpp
+++ b/src/simulation/GreensInit.cpp
@@ -121,6 +121,7 @@ void GreensInit::init(SimFramework *_sim) {
sim->console() << "# Greens normal matrix takes " << abbr_normal_bytes << " " << space_vals[norm_ind] << std::endl;
#ifdef MPI_C_FOUND
+
if (sim->getWorldSize() > 1) {
double global_shear_bytes, global_normal_bytes;
MPI_Reduce(&shear_bytes, &global_shear_bytes, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
@@ -132,5 +133,6 @@ void GreensInit::init(SimFramework *_sim) {
sim->console() << "# Global Greens shear matrix takes " << abbr_global_shear_bytes << " " << space_vals[global_shear_ind] << "." << std::endl;
sim->console() << "# Global Greens normal matrix takes " << abbr_global_normal_bytes << " " << space_vals[global_norm_ind] << "." << std::endl;
}
+
#endif
}
diff --git a/src/simulation/ProgressMonitor.cpp b/src/simulation/ProgressMonitor.cpp
index 8990160..40e113d 100644
--- a/src/simulation/ProgressMonitor.cpp
+++ b/src/simulation/ProgressMonitor.cpp
@@ -76,7 +76,7 @@ void ProgressMonitor::getStats(Simulation *sim, BlockVal &min, BlockVal &avg, Bl
// Get the minimum, maximum and sum CFF values on this node
for (lid=0; lid<sim->numLocalBlocks(); ++lid) {
gid = sim->getGlobalBID(lid);
- cur_cff = sim->getBlock(gid).getCFF();
+ cur_cff = sim->getCFF(gid);
if (cur_cff < min_val.val) {
min_val.val = cur_cff;
diff --git a/src/simulation/RunEvent.cpp b/src/simulation/RunEvent.cpp
index 670ae18..09955d6 100644
--- a/src/simulation/RunEvent.cpp
+++ b/src/simulation/RunEvent.cpp
@@ -36,10 +36,10 @@ void RunEvent::markBlocks2Fail(Simulation *sim, const FaultID &trigger_fault) {
if (sim->getBlock(gid).getFailed()) continue;
// Add this block if it has a static CFF failure
- add = sim->getBlock(gid).cffFailure();
+ add = sim->cffFailure(gid);
// Allow dynamic failure if the block is "loose" (next to a previously failed block)
- if (loose_elements.count(gid) > 0) add |= sim->getBlock(gid).dynamicFailure(trigger_fault);
+ if (loose_elements.count(gid) > 0) add |= sim->dynamicFailure(gid, trigger_fault);
if (add) {
sim->getBlock(gid).setFailed(true);
@@ -58,15 +58,16 @@ void RunEvent::processBlocksOrigFail(Simulation *sim, quakelib::ModelSweeps &swe
// For each block that fails in this sweep, calculate how much it slips
for (fit=local_failed_elements.begin(); fit!=local_failed_elements.end(); ++fit) {
if (sim->isLocalBlockID(*fit)) {
+ BlockID gid = *fit;
Block &b = sim->getBlock(*fit);
// calculate the drop in stress from the failure
- stress_drop = b.getCFF0() - b.getCFF();
+ stress_drop = sim->getCFF0(gid) - sim->getCFF(gid);
- if (!stress_drop) stress_drop = b.getStressDrop() - b.getCFF();
+ if (!stress_drop) stress_drop = sim->getStressDrop(gid) - sim->getCFF(gid);
// Slip is in m
- slip = (stress_drop/b.getSelfStresses());
+ slip = (stress_drop/sim->getSelfStresses(gid));
if (slip < 0) slip = 0;
@@ -78,10 +79,10 @@ void RunEvent::processBlocksOrigFail(Simulation *sim, quakelib::ModelSweeps &swe
b.lame_mu());
sweeps.setInitStresses(sweep_num,
b.getBlockID(),
- b.getShearStress(),
- b.getNormalStress());
+ sim->getShearStress(gid),
+ sim->getNormalStress(gid));
- b.state.slipDeficit += slip;
+ sim->setSlipDeficit(gid, sim->getSlipDeficit(gid)+slip);
}
}
}
@@ -145,17 +146,15 @@ void RunEvent::processBlocksSecondaryFailures(Simulation *sim, quakelib::ModelSw
double *x = new double[num_local_failed];
for (i=0,it=local_id_list.begin(); it!=local_id_list.end(); ++i,++it) {
- Block &blk = sim->getBlock(*it);
-
for (n=0,jt=global_id_list.begin(); jt!=global_id_list.end(); ++n,++jt) {
A[i*num_global_failed+n] = sim->getGreenShear(*it, jt->first);
if (sim->doNormalStress()) {
- A[i*num_global_failed+n] -= blk.friction()*sim->getGreenNormal(*it, jt->first);
+ A[i*num_global_failed+n] -= sim->getFriction(*it)*sim->getGreenNormal(*it, jt->first);
}
}
- b[i] = blk.getCFF()+blk.friction()*blk.getRhogd();
+ b[i] = sim->getCFF(*it)+sim->getFriction(*it)*sim->getRhogd(*it);
}
if (sim->isRootNode()) {
@@ -221,21 +220,21 @@ void RunEvent::processBlocksSecondaryFailures(Simulation *sim, quakelib::ModelSw
for (i=0,it=local_id_list.begin(); it!=local_id_list.end(); ++i,++it) {
Block &block = sim->getBlock(*it);
- double slip = x[i] - block.getSlipDeficit();
+ double slip = x[i] - sim->getSlipDeficit(*it);;
if (slip > 0) {
// Record how much the block slipped in this sweep and initial stresses
sweeps.setSlipAndArea(sweep_num,
- block.getBlockID(),
+ *it,
slip,
block.area(),
block.lame_mu());
sweeps.setInitStresses(sweep_num,
- block.getBlockID(),
- block.getShearStress(),
- block.getNormalStress());
+ *it,
+ sim->getShearStress(*it),
+ sim->getNormalStress(*it));
- block.state.slipDeficit += slip;
+ sim->setSlipDeficit(*it, sim->getSlipDeficit(*it)+slip);
}
}
@@ -290,9 +289,10 @@ void RunEvent::processStaticFailure(Simulation *sim) {
// Recalculate CFF for all blocks where slipped blocks don't contribute
for (it=sim->begin(); it!=sim->end(); ++it) {
- sim->setShearStress(it->getBlockID(), 0.0);
- sim->setNormalStress(it->getBlockID(), it->getRhogd());
- sim->setUpdateField(it->getBlockID(), (it->getFailed() ? 0 : it->state.slipDeficit));
+ BlockID gid = it->getBlockID();
+ sim->setShearStress(gid, 0.0);
+ sim->setNormalStress(gid, sim->getRhogd(gid));
+ sim->setUpdateField(gid, (it->getFailed() ? 0 : sim->getSlipDeficit(gid)));
}
// Distribute the update field values to other processors
@@ -332,8 +332,8 @@ void RunEvent::processStaticFailure(Simulation *sim) {
for (gid=0; gid<sim->numGlobalBlocks(); ++gid) {
Block &block=sim->getBlock(gid);
sim->setShearStress(gid, 0.0);
- sim->setNormalStress(gid, block.getRhogd());
- sim->setUpdateField(gid, (block.getFailed() ? 0 : block.state.slipDeficit));
+ sim->setNormalStress(gid, sim->getRhogd(gid));
+ sim->setUpdateField(gid, (block.getFailed() ? 0 : sim->getSlipDeficit(gid)));
}
sim->distributeUpdateField();
@@ -358,11 +358,10 @@ void RunEvent::processStaticFailure(Simulation *sim) {
for (fit=global_failed_elements.begin(); fit!=global_failed_elements.end(); ++fit) {
if (sim->isLocalBlockID(fit->first)) {
- Block &b = sim->getBlock(fit->first);
event_sweeps.setFinalStresses(sweep_num,
fit->first,
- b.getShearStress(),
- b.getNormalStress());
+ sim->getShearStress(fit->first),
+ sim->getNormalStress(fit->first));
}
}
@@ -439,20 +438,19 @@ void RunEvent::processAftershock(Simulation *sim) {
double element_slip = element_moment/(b.lame_mu()*b.area());
// Adjust the slip deficit appropriately
- b.state.slipDeficit += element_slip;
+ sim->setSlipDeficit(*bit, sim->getSlipDeficit(*bit)+element_slip);
// Create the sweep describing this aftershock
// Since we don't distinguish sweeps, every slip occurs in sweep 0
event_sweeps.setSlipAndArea(0, *bit, element_slip, b.area(), b.lame_mu());
- event_sweeps.setInitStresses(0, *bit, b.getShearStress(), b.getNormalStress());
+ event_sweeps.setInitStresses(0, *bit, sim->getShearStress(*bit), sim->getNormalStress(*bit));
}
// Set the update field to the slip of all blocks
for (gid=0; gid<sim->numGlobalBlocks(); ++gid) {
- Block &block=sim->getBlock(gid);
sim->setShearStress(gid, 0.0);
- sim->setNormalStress(gid, block.getRhogd());
- sim->setUpdateField(gid, block.state.slipDeficit);
+ sim->setNormalStress(gid, sim->getRhogd(gid));
+ sim->setUpdateField(gid, sim->getSlipDeficit(gid));
}
sim->distributeUpdateField();
@@ -474,8 +472,7 @@ void RunEvent::processAftershock(Simulation *sim) {
// Record final stresses on each block involved in the aftershock
for (bit=id_set.begin(); bit!=id_set.end(); ++bit) {
- Block &b=sim->getBlock(*bit);
- event_sweeps.setFinalStresses(0, *bit, b.getShearStress(), b.getNormalStress());
+ event_sweeps.setFinalStresses(0, *bit, sim->getShearStress(*bit), sim->getNormalStress(*bit));
}
// Add sweeps to list
@@ -488,7 +485,7 @@ SimRequest RunEvent::run(SimFramework *_sim) {
// Save stress information at the beginning of the event
// This is used to determine dynamic block failure
- for (lid=0; lid<sim->numLocalBlocks(); ++lid) sim->getBlock(sim->getGlobalBID(lid)).saveStresses();
+ for (lid=0; lid<sim->numLocalBlocks(); ++lid) sim->saveStresses(sim->getGlobalBID(lid));
// If there's a specific block that triggered the event, it's a static stress failure type event
if (sim->getCurrentEvent().getEventTrigger() != UNDEFINED_ELEMENT_ID) {
diff --git a/src/simulation/SanityCheck.cpp b/src/simulation/SanityCheck.cpp
index dcfa722..013bcce 100644
--- a/src/simulation/SanityCheck.cpp
+++ b/src/simulation/SanityCheck.cpp
@@ -33,7 +33,7 @@ bool SanityCheck::assertCFFValueCorrectness(Simulation *sim) {
for (lid=0; lid<sim->numLocalBlocks(); ++lid) {
gid = sim->getGlobalBID(lid);
- double val = sim->getBlock(gid).getCFF();
+ double val = sim->getCFF(gid);
if (std::isnan(val) || fabs(val) > 1e20) {
failed_cffs.push_back(std::make_pair(gid, val));
diff --git a/src/simulation/UpdateBlockStress.cpp b/src/simulation/UpdateBlockStress.cpp
index 3e11573..6d54985 100644
--- a/src/simulation/UpdateBlockStress.cpp
+++ b/src/simulation/UpdateBlockStress.cpp
@@ -36,6 +36,12 @@ void UpdateBlockStress::init(SimFramework *_sim) {
for (lid=0; lid<sim->numLocalBlocks(); ++lid) {
gid = sim->getGlobalBID(lid);
+ double rho = 5.515e3; // density of rock in kg m^-3
+ double g = 9.81; // force of gravity in m s^-2
+ double depth = -sim->getBlock(gid).center()[2]; // depth of block center in m
+
+ sim->setRhogd(gid, rho*g*depth); // kg m^-3 * m s^-2 * m = kg m^-1 * s^-2 = Pa
+
// Set stresses to their specified initial values
sim->setShearStress(gid, sim->getBlock(gid).getInitShearStress());
sim->setNormalStress(gid, sim->getBlock(gid).getInitNormalStress());
@@ -48,10 +54,10 @@ void UpdateBlockStress::init(SimFramework *_sim) {
stress_drop += (nt->slip_rate()/norm_velocity)*sim->getGreenShear(gid, nt->getBlockID());
}
- sim->getBlock(gid).setStressDrop(sim->getBlock(gid).max_slip()*stress_drop);
+ sim->setStressDrop(gid, sim->getBlock(gid).max_slip()*stress_drop);
// noise in the range [1-a, 1+a)
- sim->getBlock(gid).state.slipDeficit = 0;
+ sim->setSlipDeficit(gid, 0);
if (sim->isLocalBlockID(gid)) {
sim->decompressNormalRow(gid);
@@ -107,8 +113,10 @@ SimRequest UpdateBlockStress::run(SimFramework *_sim) {
sim->incrementYear(next_event_global.val);
for (lid=0; lid<sim->numLocalBlocks(); ++lid) {
- Block &local_block = sim->getBlock(sim->getGlobalBID(lid));
- local_block.state.slipDeficit -= local_block.slip_rate()*convert.year2sec(next_event_global.val)*(1.0-local_block.aseismic());
+ BlockID gid = sim->getGlobalBID(lid);
+ Block &local_block = sim->getBlock(gid);
+ double cur_slip_deficit = sim->getSlipDeficit(gid);
+ sim->setSlipDeficit(gid, cur_slip_deficit-local_block.slip_rate()*convert.year2sec(next_event_global.val)*(1.0-local_block.aseismic()));
}
// Recalculate the stress on all blocks with the new slip deficits
@@ -166,7 +174,8 @@ void UpdateBlockStress::nextStaticFailure(BlockVal &next_static_fail) {
if (sim->doNormalStress()) {
for (it=sim->begin(); it!=sim->end(); ++it) {
- sim->setUpdateField(it->getBlockID(), -it->friction()*it->slip_rate());
+ BlockID gid = it->getBlockID();
+ sim->setUpdateField(gid, -sim->getFriction(gid)*it->slip_rate());
}
sim->matrixVectorMultiplyAccum(tmpBuffer,
@@ -189,11 +198,11 @@ void UpdateBlockStress::nextStaticFailure(BlockVal &next_static_fail) {
if (block.aseismic() > 0) {
double A, B, K;
A = -tmpBuffer[gid];
- B = -block.aseismic()*block.getSelfStresses()/block.getRecurrence();
- K = -log(A+B*block.getCFF())/B;
+ B = -block.aseismic()*sim->getSelfStresses(gid)/sim->getRecurrence(gid);
+ K = -log(A+B*sim->getCFF(gid))/B;
ts = K + log(A)/B;
} else {
- ts = convert.sec2year(block.getCFF()/tmpBuffer[gid]);
+ ts = convert.sec2year(sim->getCFF(gid)/tmpBuffer[gid]);
}
// Blocks with negative timesteps are skipped. These effectively mean the block
@@ -223,9 +232,10 @@ void UpdateBlockStress::stressRecompute(void) {
// Reset the shear and normal stress, and set the update field to be the slip deficit
for (it=sim->begin(); it!=sim->end(); ++it) {
- sim->setShearStress(it->getBlockID(), 0.0);
- sim->setNormalStress(it->getBlockID(), it->getRhogd());
- sim->setUpdateField(it->getBlockID(), it->state.slipDeficit);
+ BlockID gid = it->getBlockID();
+ sim->setShearStress(gid, 0.0);
+ sim->setNormalStress(gid, sim->getRhogd(gid));
+ sim->setUpdateField(gid, sim->getSlipDeficit(gid));
}
// Distribute the new update field over all nodes
More information about the CIG-COMMITS
mailing list