Muskingum musle

crhmcode/Makefile

@@ -13,7 +13,7 @@ CFLAGS = -Wno-switch -Wno-c++11-extensions -Wno-logical-op-parentheses \
 
 CFLAGS = -I. -I$(BOOST) -I$(SDIR)/core -I$(SDIR)/macro \
 -I$(SDIR)/modules/waterquality -I$(SDIR)/modules/newmodules \
--I$(SDIR)/modules/coremodules -I$(SDIR)/gcc -std=c++14 -static -static-libgcc -static-libstdc++
+-I$(SDIR)/modules/coremodules -I$(SDIR)/gcc -std=c++14 -static-libgcc -static-libstdc++ -lfmt
 
 CORE_OBJS = GlobalDll.o StandardConverterUtility.o
 
@@ -30,7 +30,7 @@ Classramp.o Classrandom.o Classrefwind.o ClassReplace.o Classrh.o \
 ClassRH_WtoI.o ClassSim.o Classsin.o ClassSmear.o Classsquare.o Classsub.o \
 ClasssubV.o Classtime.o ClassTimeshift.o
 
-ClassModule_OBJS = Administer.o ClassClark.o ClassModule.o ClassMuskingum.o Myparser.o
+ClassModule_OBJS = Administer.o ClassClark.o ClassModule.o ClassMuskingum.o ClassMuskingum2.o Myparser.o
 
 Common_OBJS = Common.o snowcover.o
 
@@ -61,7 +61,7 @@ ClassGrow_Crop.o ClassHMSA.o ClassHtobs.o ClassIceBulb.o ClassICEflow.o \
 ClassIntcp.o Classinterception.o ClassK_Estimate.o ClassKevin.o Classlake.o \
 ClassLongVt.o ClassMeltRunoff_Kstorage.o ClassMeltRunoff_Lag.o ClassMod_Exec.o \
 ClassNeedle.o Classnetall.o ClassNetroute.o ClassNetroute_D.o ClassNetroute_M.o \
-ClassNetroute_M_D.o ClassNO_pbsm.o ClassNOP.o ClassObs.o ClassObstoPar.o \
+ClassNetroute_M_D.o ClassNetroute_M2_D.o ClassNO_pbsm.o ClassNOP.o ClassObs.o ClassObstoPar.o \
 Classpbsm.o Classpbsm_M.o ClasspbsmSnobal.o ClassPrairieInfil.o ClassPSPnew.o \
 Classqdrift.o Classqmelt.o Classquinton.o ClassREWroute.o ClassREWroute_stream.o \
 ClassREWroute2.o Classsbsm.o ClassSetSoil.o Classshared.o ClassShutWall.o \
@@ -91,7 +91,7 @@ OBJ = $(patsubst %,$(ODIR)/%,$(ALL_OBS))
 crhm: $(OBJ)
 	$(CC) -o $@ $^ $(CFLAGS) $(DFLAGS) $(LIBS)
 
-$(ODIR)/main.o: $(SDIR)/gcc/main.cpp | $(ODIR)
+$(ODIR)/main.o: $(SDIR)/main.cpp | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
 $(ODIR)/%.o: $(SDIR)/core/%.cpp $(SDIR)/core/%.h | $(ODIR)
@@ -103,7 +103,7 @@ $(ODIR)/%.o: $(SDIR)/core/ClassModule/%.cpp $(SDIR)/core/ClassModule/%.h | $(ODI
 $(ODIR)/%.o: $(SDIR)/core/ClassCRHM/%.cpp $(SDIR)/core/ClassCRHM/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
-$(ODIR)/%.o: $(SDIR)/core/ClassCRHM/Filters/%.cpp  $(SDIR)/core/ClassCRHM/Filters/%.h | $(ODIR)
+$(ODIR)/%.o: $(SDIR)/core/filter/%.cpp  $(SDIR)/core/filter/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
 $(ODIR)/%.o: $(SDIR)/core/TStringList/%.cpp $(SDIR)/core/TStringList/%.h | $(ODIR)
@@ -112,10 +112,10 @@ $(ODIR)/%.o: $(SDIR)/core/TStringList/%.cpp $(SDIR)/core/TStringList/%.h | $(ODI
 $(ODIR)/%.o: $(SDIR)/core/Common/%.cpp $(SDIR)/core/Common/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
-$(ODIR)/%.o: $(SDIR)/core/CRHM_parse/%.cpp $(SDIR)/core/CRHM_parse/%.h | $(ODIR)
+$(ODIR)/%.o: $(SDIR)/core/def/%.cpp $(SDIR)/core/def/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
-$(ODIR)/%.o: $(SDIR)/core/CRHM_parse/exec/%.cpp $(SDIR)/core/CRHM_parse/exec/%.h | $(ODIR)
+$(ODIR)/%.o: $(SDIR)/core/exec/%.cpp $(SDIR)/core/exec/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
 $(ODIR)/%.o: $(SDIR)/core/CRHMmain/%.cpp $(SDIR)/core/CRHMmain/%.h | $(ODIR)
@@ -136,7 +136,7 @@ $(ODIR)/%.o: $(SDIR)/modules/waterquality/%.cpp $(SDIR)/modules/waterquality/%.h
 $(ODIR)/%.o: $(SDIR)/modules/newmodules/%.cpp $(SDIR)/modules/newmodules/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
-$(ODIR)/%.o: $(SDIR)/modules/coremodules/%.cpp $(SDIR)/modules/coremodules/%.h | $(ODIR)
+$(ODIR)/%.o: $(SDIR)/modules/%.cpp $(SDIR)/modules/%.h | $(ODIR)
 	$(CC) -c -o $@ $< $(CFLAGS) $(DFLAGS)
 
 clean:

crhmcode/src/core/CMakeLists.txt

@@ -27,6 +27,7 @@ add_library(core OBJECT
     ClassModel.cpp
     ClassModule.cpp
     ClassMuskingum.cpp
+    ClassMuskingum2.cpp
     ClassPar.cpp
     ClassVar.cpp
     Common.cpp

crhmcode/src/core/ClassModule.h

@@ -17,6 +17,7 @@
 #include "Myparser.h"
 #include "ClassClark.h"
 #include "ClassMuskingum.h"
+#include "ClassMuskingum2.h"
 
 
 

crhmcode/src/core/ClassMuskingum2.cpp

@@ -0,0 +1,109 @@
+#include "ClassMuskingum2.h"
+#include <algorithm>    // std::max
+#include <exception>    // std::exception
+
+
+//---------------------------------------------------------------------------
+ClassMuskingum2::ClassMuskingum2(const double* inVar, double* outVar, 
+                                const double* k, const double* X_M, const long nhru)
+	: inVar(inVar), outVar(outVar), nhru(nhru) {
+
+	// !!! UNITS !!!
+
+	// kstorage (days)
+	// lag (hours)
+
+	c0 = new double[nhru];
+	c1 = new double[nhru];
+	c2 = new double[nhru];
+
+	/*
+	* If on every cycle the water moves one step then the number of elements determines how many steps
+	* are required to get to the end of the stream.
+	*/
+
+	buff_sz = new long[nhru];
+	K_sub = new double[nhru];
+	x_atten = new double[nhru];
+
+
+	for (long hh = 0; hh < nhru; hh++) {
+
+// Add 1 to buff_sz because the first buffer location just holds the inflow values
+// actual outflow values start from the second buffer location
+        buff_sz[hh] = 1 + std::max(1l, lround(k[hh]/Global::Interval));
+        K_sub[hh] = k[hh] / ( buff_sz[hh] - 1);  // The 1 additional should not be included here
+        x_atten[hh] = X_M[hh];
+
+
+		c0[hh] = (Global::Interval - 2.0 * K_sub[hh] * X_M[hh]) /
+			(2.0 * K_sub[hh] * (1.0 - X_M[hh]) + Global::Interval);  // units of Global::Interval (days)
+
+		c1[hh] = (Global::Interval + 2.0 * K_sub[hh] * X_M[hh]) /
+			(2.0 * K_sub[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days)
+
+		c2[hh] = (2.0 * K_sub[hh] * (1.0 - X_M[hh]) - Global::Interval) /
+			(2.0 * K_sub[hh] * (1.0 - X_M[hh]) + Global::Interval); // units of kstorage (days)
+
+	}
+
+	buff_q = new double* [nhru];   // create mem storage for segment discharge
+
+	for (long hh = 0; hh < nhru; hh++) {
+		buff_q[hh] = new double[buff_sz[hh]];
+		for (long jj = 0; jj < buff_sz[hh]; jj++)
+			buff_q[hh][jj] = 0.0;
+	}
+}
+
+ClassMuskingum2::~ClassMuskingum2() {
+	delete[] c0;
+	delete[] c1;
+	delete[] c2;
+	delete[] buff_sz;
+	delete[] K_sub;
+	delete[] x_atten;
+
+	for (long hh = 0; hh < nhru; hh++)
+		delete[] buff_q[hh];
+	delete[] buff_q;
+}
+
+void ClassMuskingum2::ChangeLag(const double* newlag, const long hh)
+{
+    // This is only ever used by ClassMeltRunoff_Lag
+    throw std::runtime_error("Not Implemented");
+}
+
+void ClassMuskingum2::DoMuskingum() {
+
+	for (long hh = 0; hh < nhru; hh++) {
+        DoMuskingum(hh);
+	}
+}
+
+void ClassMuskingum2::DoMuskingum(const long hh) {
+
+    double q_prev = buff_q[hh][0];
+    buff_q[hh][0] = inVar[hh];
+    for (int i=1; i<buff_sz[hh]; i++) {
+        const double q_tmp = buff_q[hh][i];   
+    // Qout1 = C0 * Qin1 + C1 * Qin0 + C2 * Qout0 
+        buff_q[hh][i] = c0[hh] * buff_q[hh][i-1] + c1[hh] * q_prev + c2[hh] * buff_q[hh][i];
+        q_prev = q_tmp;
+    }
+    outVar[hh] = buff_q[hh][buff_sz[hh]-1];
+}
+
+double ClassMuskingum2::Left(int hh) {
+
+	double Sstorage = 0;
+
+	for (int ii = 1; ii < buff_sz[hh]; ++ii) {
+        Sstorage += K_sub[hh] * ( x_atten[hh] * buff_q[hh][ii-1] + (1-x_atten[hh])*buff_q[hh][ii] );
+    }
+
+	return Sstorage;
+}
+
+//---------------------------------------------------------------------------

crhmcode/src/core/ClassMuskingum2.h

@@ -0,0 +1,40 @@
+#pragma once
+
+#ifndef CLASSMUSKINGUM2
+#define CLASSMUSKINGUM2
+
+#include "ClassModule.h"
+
+
+class   ClassMuskingum2 {
+
+public:
+	ClassMuskingum2(const double* inVar, double* outVar, 
+                    const double* k, const double* X_M, const long nhru);
+	~ClassMuskingum2();
+	void DoMuskingum();
+	void DoMuskingum(const long hh);
+	void ChangeLag(const double* newlag, const long hh);
+	double Left(int hh);
+
+
+	double* c0; // storage constant from K
+	double* c1; // storage constant from K
+	double* c2; // storage constant from K
+	double prevdate{ 0.0 };
+
+private:
+	const  double* inVar;
+	double* outVar;
+
+	double** buff_q;
+
+	double* K_sub; //
+	double* x_atten; //
+
+	long nhru;
+	long* buff_sz; // number of segments (buckets) 
+
+};
+
+#endif // !ClassMuskingum2

crhmcode/src/modules/CMakeLists.txt

@@ -55,6 +55,7 @@ add_library(modules OBJECT
     Classnetall.cpp
     ClassNetroute_D.cpp
     ClassNetroute_M_D.cpp
+    ClassNetroute_M2_D.cpp
     ClassNetroute_M.cpp
     ClassNetroute.cpp
     ClassNO_pbsm.cpp