From 85b900363d58b42aa0d0dc4d5e3a2d02bf792448 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 16 Jun 2020 12:43:02 +0200 Subject: [PATCH] =?UTF-8?q?Block=20decomposition:=20for=20=E2=80=9Csolve?= =?UTF-8?q?=20two=20boundaries=E2=80=9D=20block,=20no=20longer=20compute?= =?UTF-8?q?=20=E2=80=9Cb=E2=80=9D=20(-residuals+g1*y)=20in=20the=20dynamic?= =?UTF-8?q?=20function?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/DynamicModel.cc | 31 +++---------------------------- 1 file changed, 3 insertions(+), 28 deletions(-) diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index b330bcd8..a03ceb03 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -244,8 +244,6 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const ofstream output; output.open(filename, ios::out | ios::binary); - vector Uf(block_size); - output << "%" << endl << "% " << filename << " : Computes dynamic version of one block" << endl << "%" << endl @@ -261,7 +259,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const || simulation_type == BlockSimulationType::solveForwardSimple) output << "function [residual, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, it_, jacobian_eval)" << endl; else - output << "function [residual, T, g1, b, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, periods, jacobian_eval, y_kmin, y_size)" << endl; + output << "function [residual, T, g1, varargout] = dynamic_" << blk+1 << "(y, x, params, steady_state, T, periods, jacobian_eval, y_kmin, y_size)" << endl; output << " % ////////////////////////////////////////////////////////////////////////" << endl << " % //" << string(" Block ").substr(static_cast(log10(blk + 1))) << blk+1 @@ -310,8 +308,7 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const if (simulation_type == BlockSimulationType::solveTwoBoundariesComplete || simulation_type == BlockSimulationType::solveTwoBoundariesSimple) - output << " b = zeros(periods*y_size,1);" << endl - << " for it_ = y_kmin+1:(periods+y_kmin)" << endl + output << " for it_ = y_kmin+1:(periods+y_kmin)" << endl << " Per_y_=it_*y_size;" << endl << " Per_J_=(it_-y_kmin-1)*y_size;" << endl << " Per_K_=(it_-1)*y_size;" << endl; @@ -381,7 +378,6 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const case BlockSimulationType::solveTwoBoundariesSimple: if (eq < block_recursive_size) goto evaluation; - Uf[eq] << "b(" << eq+1-block_recursive_size << "+Per_J_) = -residual(" << eq+1-block_recursive_size << ", it_)"; output << indent_prefix << "residual(" << eq+1-block_recursive_size << ", it_) = ("; goto end; default: @@ -466,26 +462,8 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const for (const auto &[indices, d] : blocks_derivatives[blk]) { auto [eq, var, lag] = indices; - int varr = getBlockVariableID(blk, var); if (eq >= block_recursive_size && var >= block_recursive_size) { - if (lag == 0) - Uf[eq] << "+g1(" << eq+1-block_recursive_size - << "+Per_J_, " << var+1-block_recursive_size - << "+Per_K_)*y(it_, " << varr+1 << ")"; - else if (lag == 1) - Uf[eq] << "+g1(" << eq+1-block_recursive_size - << "+Per_J_, " << var+1-block_recursive_size - << "+Per_y_)*y(it_+1, " << varr+1 << ")"; - else if (lag > 0) - Uf[eq] << "+g1(" << eq+1-block_recursive_size - << "+Per_J_, " << var+1-block_recursive_size - << "+y_size*(it_+" << lag-1 << "))*y(it_+" << lag << ", " << varr+1 << ")"; - else - Uf[eq] << "+g1(" << eq+1-block_recursive_size - << "+Per_J_, " << var+1-block_recursive_size - << "+y_size*(it_" << lag-1 << "))*y(it_" << lag << ", " << varr+1 << ")"; - output << " "; if (lag == 0) output << " g1(" << eq+1-block_recursive_size << "+Per_J_, " @@ -503,9 +481,6 @@ DynamicModel::writeDynamicPerBlockMFiles(const string &basename) const output << ";" << endl; } } - for (int eq = 0; eq < block_size; eq++) - if (eq >= block_recursive_size) - output << " " << Uf[eq].str() << ";" << endl; output << " end" << endl << " end" << endl; @@ -1392,7 +1367,7 @@ DynamicModel::writeDynamicBlockMFile(const string &basename) const break; case BlockSimulationType::solveTwoBoundariesComplete: case BlockSimulationType::solveTwoBoundariesSimple: - mDynamicModelFile << " [r, T, dr(" << blk + 1 << ").g1, b, dr(" << blk + 1 << ").g1_x, dr(" << blk + 1 << ").g1_xd, dr(" << blk + 1 << ").g1_o]=" << basename << ".block.dynamic_" << blk + 1 << "(y, x, params, steady_state, T, it_-" << max_lag << ", true, " << max_lag << ", " << block_recursive_size << ");" << endl + mDynamicModelFile << " [r, T, dr(" << blk + 1 << ").g1, dr(" << blk + 1 << ").g1_x, dr(" << blk + 1 << ").g1_xd, dr(" << blk + 1 << ").g1_o]=" << basename << ".block.dynamic_" << blk + 1 << "(y, x, params, steady_state, T, it_-" << max_lag << ", true, " << max_lag << ", " << block_recursive_size << ");" << endl << " residual(y_index_eq)=r(:,M_.maximum_lag+1);" << endl; break; default: