Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MOD files containing two SOLVEs don't compile. #1501

Open
1uc opened this issue Oct 7, 2024 · 1 comment
Open

MOD files containing two SOLVEs don't compile. #1501

1uc opened this issue Oct 7, 2024 · 1 comment

Comments

@1uc
Copy link
Collaborator

1uc commented Oct 7, 2024

For example the following:

$ cat two_blocks.mod
NEURON {
    SUFFIX two_blocks
    RANGE A, B, X, Y
    NONSPECIFIC_CURRENT il
}

STATE {
    X
    Y
    A
    B
}

ASSIGNED {
    il
    i
}

BREAKPOINT {
    SOLVE state1 METHOD sparse
    SOLVE state2 METHOD sparse
    il = i
}

INITIAL {
    A = 1.0
    B = 3.0
    X = 0.0
    Y = 1.0
}

KINETIC state1 {
    ~ X <-> Y (0.4, 0.5)
}

KINETIC state2 {
    ~ A <-> B (0.6, 0.7)
    i = (f_flux-b_flux)
}

generates the following C++ code:

            Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
            double* nmodl_eigen_x = nmodl_eigen_xm.data();
            nmodl_eigen_x[static_cast<int>(0)] = inst->X[id];
            nmodl_eigen_x[static_cast<int>(1)] = inst->Y[id];
            // call newton solver
            functor_two_blocks_0 newton_functor(nt, inst, id, pnodecount, v, indexes, data, thread);
            newton_functor.initialize();
            int newton_iterations = nmodl::newton::newton_solver(nmodl_eigen_xm, newton_functor);
            if (newton_iterations < 0) assert(false && "Newton solver did not converge!");
            inst->X[id] = nmodl_eigen_x[static_cast<int>(0)];
            inst->Y[id] = nmodl_eigen_x[static_cast<int>(1)];
            newton_functor.initialize(); // TODO mimic calling F again.
            newton_functor.finalize();

            
            Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
            double* nmodl_eigen_x = nmodl_eigen_xm.data();
            nmodl_eigen_x[static_cast<int>(0)] = inst->A[id];
            nmodl_eigen_x[static_cast<int>(1)] = inst->B[id];
            // call newton solver
            functor_two_blocks_1 newton_functor(nt, inst, id, pnodecount, v, indexes, data, thread);
            newton_functor.initialize();
            int newton_iterations = nmodl::newton::newton_solver(nmodl_eigen_xm, newton_functor);
            if (newton_iterations < 0) assert(false && "Newton solver did not converge!");
            inst->A[id] = nmodl_eigen_x[static_cast<int>(0)];
            inst->B[id] = nmodl_eigen_x[static_cast<int>(1)];
            newton_functor.initialize(); // TODO mimic calling F again.
            newton_functor.finalize();

and fails with:

$ nrnivmodl -coreneuron
x86_64/corenrn/mod2c/two_blocks.cpp:786:41: error: redeclaration of ‘Eigen::Matrix<double, 2, 1> nmodl_eigen_xm’
  786 |             Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
      |                                         ^~~~~~~~~~~~~~
x86_64/corenrn/mod2c/two_blocks.cpp:771:41: note: ‘Eigen::Matrix<double, 2, 1> nmodl_eigen_xm’ previously declared here
  771 |             Eigen::Matrix<double, 2, 1> nmodl_eigen_xm;
      |                                         ^~~~~~~~~~~~~~

NOCMODL seems to generate code that compiles and looks sane. Therefore, it's not obvious the MOD file is invalid.

@1uc
Copy link
Collaborator Author

1uc commented Oct 7, 2024

Summary of discussion with @nrnhines:

  1. Carefully check if NOCMODL supports this and check ModelDB CI for usage examples.
  2. Also consider recursive use, e.g. solving a LINEAR block inside a KINETIC block.
  3. Note that STATE variables are declared globally, i.e. there's no notion of the STATEs of block state1. Some solvers might require knowing the number of states.

Generally, this seems to be a problem that can be revisited when there's an instance of this problem.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant