partial.cpp

Go to the documentation of this file.

#include <../../nmodlconf.h>
/* /local/src/master/nrn/src/nmodl/partial.c,v 4.1 1997/08/30 20:45:33 hines Exp */
/*
partial.c,v
 * Revision 4.1  1997/08/30  20:45:33  hines
 * cvs problem with branches. Latest nmodl stuff should now be a top level
 *
 * Revision 4.0.1.1  1997/08/08  17:24:00  hines
 * nocmodl version 4.0.1
 *
 * Revision 4.0  1997/08/08  17:06:28  hines
 * proper nocmodl version number
 *
 * Revision 1.1.1.1  1994/10/12  17:21:37  hines
 * NEURON 3.0 distribution
 *
 * Revision 9.76  90/12/07  09:27:23  hines
 * new list structure that uses unions instead of void *element
 *
 * Revision 9.67  90/11/27  12:10:07  hines
 * bug in bnd introduced by change to partial
 *
 * Revision 9.66  90/11/27  10:47:26  hines
 * allow multiple partial equations within a partial block
 *
 * Revision 9.58  90/11/20  17:24:15  hines
 * CONSTANT changed to PARAMETER
 * CONSTANT now refers to variables that don't get put in .var file
 *
 * Revision 9.20  90/08/15  15:16:16  hines
 * need to pass error value from crank back through the
 * block called by SOLVE.
 *
 * Revision 8.1  90/01/16  11:06:14  mlh
 * error checking and cleanup after error and call to abort_run()
 *
 * Revision 8.0  89/09/22  17:26:54  nfh
 * Freezing
 *
 * Revision 7.0  89/08/30  13:32:29  nfh
 * Rev 7 is now Experimental; Rev 6 is Testing
 *
 * Revision 6.0  89/08/14  16:27:10  nfh
 * Rev 6.0 is latest of 4.x; now the Experimental version
 *
 * Revision 4.1  89/07/27  13:27:30  mlh
 * crank calling sequence sends indepsym as a sentinal value
 * ugh!
 *
 * Revision 4.0  89/07/24  17:03:40  nfh
 * Freezing rev 3.  Rev 4 is now Experimental
 *
 * Revision 3.3  89/07/18  11:55:16  mlh
 * first_time removed and MODEL_LEVEL used for declaration precedence
 *
 * Revision 1.4  89/07/18  11:22:16  mlh
 * eliminate first_time, etc.
 *
 * Revision 1.3  89/07/11  16:51:58  mlh
 * remove p array from calling sequence of crank
 *
 * Revision 1.2  89/07/07  08:55:01  mlh
 * START field in INDEPENDENT block
 * SWEEP keyword replaces SCOP in INDEPENDENT block
 * FIRST and LAST keywords in PARTIAL block for boundaries rplace index
 * y'0 etc. allowed for constants
 *
 * Revision 1.1  89/07/06  14:50:25  mlh
 * Initial revision
 *
*/
 
#include "modl.h"
#include "parse1.hpp"
#include "symbol.h"
 
extern Symbol* indepsym;
 
/* for partial block communication with solve */
static Symbol* pv[4]; /* DV, F, V, G */
static Item* partialcolon;
static List* parinfo; /* DV, F, V, G, ~ with ident */
 
void solv_partial(Item* qsol, Symbol* fun) {
    int i, ident;
    Item* q;
    Symbol* dspace;
 
    if ((dspace = lookup("delta_x")) == SYM0) {
        dspace = install("delta_x", NAME);
        parminstall(dspace, "1", "", "");
    }
    Sprintf(buf, "%s();\n", fun->name);
    replacstr(qsol, buf);
    save_dt(partialcolon);
    ITERATE(q, parinfo) {
        for (i = 0; i < 4; i++) {
            pv[i] = SYM(q);
            q = q->next;
        }
        partialcolon = ITM(q);
        ident = q->itemtype;
 
        Sprintf(buf,
                "if (error=crank(%d, %s, %s, %s, delta_%s, %s, %s, _pbound%d, &_parwork%d))return "
                "error;\n",
                pv[2]->araydim,
                pv[2]->name,
                pv[1]->name,
                pv[3]->name,
                indepsym->name,
                dspace->name,
                indepsym->name,
                ident,
                ident);
        replacstr(partialcolon, buf);
        Sprintf(buf,
                "static double *_pbound%d[4], _bnd%d_0, _bnd%d_N, *_parwork%d;\n",
                ident,
                ident,
                ident,
                ident);
        Linsertstr(procfunc, buf);
    }
}
 
void partial_eqn(Item* q2, Item* q4, Item* q8, Item* q11) /*V' F V G*/
{
    int i, dim;
    Item* q;
    static int ident = 1;
 
    if (!parinfo) {
        parinfo = newlist();
    }
    pv[0] = SYM(q2);
    pv[1] = SYM(q4);
    pv[2] = SYM(q8);
    pv[3] = SYM(q11);
    ITERATE(q, parinfo) {
        if (SYM(q) == SYM(q2)) {
            diag("Two partial equations for same state: ", SYM(q8)->name);
        }
        q = q->next->next->next->next;
    }
    Lappendsym(parinfo, SYM(q2));
    Lappendsym(parinfo, SYM(q4));
    Lappendsym(parinfo, SYM(q8));
    Lappendsym(parinfo, SYM(q11));
    q = lappendsym(parinfo, SYM0);
    ITM(q) = q2->prev;
    q->itemtype = ident++;
 
    dim = pv[2]->araydim;
    if (!(pv[2]->subtype & STAT)) {
        diag(pv[2]->name, "is not a STATE");
    }
    if (strcmp(pv[2]->name, pv[0]->name + 1)) {
        diag(pv[2]->name, "' must be the time derivative");
    }
    for (i = 0; i < 4; i++) {
        if (!(pv[i]->subtype & ARRAY) || (pv[i]->araydim != dim)) {
            diag(pv[i]->name, "dimension differs from STATE var");
        }
    }
    replacstr(q2->prev, "\n/*where partial call goes*/\n");
    partialcolon = q2->prev;
    deltokens(q2, q11);
}
 
static List* bndinfo; /* symbol, expression list with itemtype 0-3 */
static List* bnd[4];  /* for each condition consisting of
              symbol, list of tokens for expression
          in the order DEL y[0]   DEL  y[N]   y[0]   y[N]*/
 
void massagepartial(Item* q1, Item* q2, Item* q3, Item* q6) /*PARTIAL NAME stmtlist '}'*/
{
    int i, ident;
    Item *q, *qb;
    Symbol* s;
 
    static int seepartial = 0;
    if (seepartial++) {
        diag("Only one partial block currently allowed", (char*) 0);
    }
    if (!pv[0]) {
        diag("within the PARTIAL block must occur at least one equation with the syntax ---\n",
             "~ V' = F*DEL2(V) + G\n");
    }
 
    if (!bndinfo) {
        bndinfo = newlist();
    }
    ITERATE(q, parinfo) {
        for (i = 0; i < 4; i++) {
            pv[i] = SYM(q);
            q = q->next;
            bnd[i] = (List*) 0;
        }
        partialcolon = ITM(q);
        ident = q->itemtype;
        /*look through bndinfo list to reconstruct bnd*/
        ITERATE(qb, bndinfo) {
            s = SYM(qb);
            qb = qb->next;
            if (s == pv[2]) {
                if (bnd[qb->itemtype]) {
                    diag("Duplicate boundary condition for ", s->name);
                }
                bnd[qb->itemtype] = LST(qb);
                qb->itemtype = -1; /* mark used */
            }
        }
 
        /* ensure that there is one condition on each side */
        for (i = 0; i < 2; i++) {
            if (bnd[i]) {         /* Neumann */
                if (bnd[i + 2]) { /* and also dirichlet */
                    diag("Neumann and Dirichlet conditions", "specified on same side");
                }
            } else {               /* Neumann not specified */
                if (!bnd[i + 2]) { /* neither is dirichlet */
                    /* then default Neumann = 0 */
                    bnd[i] = newlist();
                    Lappendstr(bnd[i], "0");
                }
            }
        }
        /* place the boundary conditions before the call to crank */
        /* and set up the pointers */
        for (i = 0; i < 4; i++) {
            if (bnd[i]) {
                if (i % 2) {
                    Sprintf(buf, "_bnd%d_N = ", ident);
                    Insertstr(partialcolon, buf);
                    Sprintf(buf, "_pbound%d[%d] = &_bnd%d_N;\n", ident, i, ident);
                } else {
                    Sprintf(buf, "_bnd%d_0 = ", ident);
                    Insertstr(partialcolon, buf);
                    Sprintf(buf, "_pbound%d[%d] = &_bnd%d_0;\n", ident, i, ident);
                }
                Lappendstr(initlist, buf);
                move(bnd[i]->next, bnd[i]->prev, partialcolon);
                Insertstr(partialcolon, ";\n");
            }
        }
    }
    /* check that all boundary conditions have been used */
    ITERATE(qb, bndinfo) {
        s = SYM(qb);
        qb = qb->next;
        if (qb->itemtype != -1) {
            diag("Boundary condition: No partial equation for ", s->name);
        }
    }
    replacstr(q1, "int");
    Insertstr(q3, "()\n");
    SYM(q2)->subtype |= PARF;
    SYM(q2)->u.i = 0; /* not used: but see listnum in solvhandler */
    Insertstr(q6, "return 0;\n");
    movelist(q1, q6, procfunc);
}
 
/* ~ optionalDEL var[index] = expr */
/* type 0 Dirichlet (no DEL), type 1 Neumann (with DEL) */
void partial_bndry(int type, Item* qvar, Item* qfirstlast, Item* qexpr, Item* qlast) {
    int indx;
    Item* q;
    List* l;
 
    if (!bndinfo) {
        bndinfo = newlist();
    }
    if (SYM(qfirstlast)->type == FIRST) {
        indx = type;
    } else {
        indx = type + 1;
    }
    Lappendsym(bndinfo, SYM(qvar));
    q = lappendsym(bndinfo, SYM0);
    q->itemtype = indx;
    l = newlist();
    LST(q) = l;
 
    if (indx < 2) { /* ~ DEL y [N] = */
        deltokens(qvar->prev->prev, qexpr->prev);
    } else { /* ~ ... = */
        deltokens(qvar->prev, qexpr->prev);
    }
    movelist(qexpr, qlast, l);
}