/[ascend]/trunk/base/generic/integrator/ida.c
ViewVC logotype

Diff of /trunk/base/generic/integrator/ida.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1394 by jpye, Sat Apr 7 14:43:31 2007 UTC revision 1395 by jpye, Sat Apr 21 15:11:45 2007 UTC
# Line 161  typedef struct IntegratorIdaPrecDJStruct Line 161  typedef struct IntegratorIdaPrecDJStruct
161      N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */      N_Vector PIii; /**< diagonal elements of the inversed Jacobi preconditioner */
162  } IntegratorIdaPrecDataJacobi;  } IntegratorIdaPrecDataJacobi;
163    
164    typedef struct IntegratorIdaPrecDJFStruct{
165        linsolqr_system_t L;
166    } IntegratorIdaPrecDataJacobian;
167    
168  /**  /**
169      Hold all the function pointers associated with a particular preconditioner      Hold all the function pointers associated with a particular preconditioner
170      We don't need to store the 'pfree' function here as it is allocated to the enginedata struct      We don't need to store the 'pfree' function here as it is allocated to the enginedata struct
# Line 222  void integrator_ida_write_stats(Integrat Line 226  void integrator_ida_write_stats(Integrat
226  void integrator_ida_write_incidence(IntegratorSystem *sys);  void integrator_ida_write_incidence(IntegratorSystem *sys);
227    
228  /*------  /*------
229      Full jacobian preconditioner -- experimental
230    */
231    
232    int integrator_ida_psetup_jacobian(realtype tt,
233             N_Vector yy, N_Vector yp, N_Vector rr,
234             realtype c_j, void *prec_data,
235             N_Vector tmp1, N_Vector tmp2,
236             N_Vector tmp3
237    );
238    
239    int integrator_ida_psolve_jacobian(realtype tt,
240             N_Vector yy, N_Vector yp, N_Vector rr,
241             N_Vector rvec, N_Vector zvec,
242             realtype c_j, realtype delta, void *prec_data,
243             N_Vector tmp
244    );
245    
246    void integrator_ida_pcreate_jacobian(IntegratorSystem *sys);
247    
248    void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata);
249    
250    static const IntegratorIdaPrec prec_jacobian = {
251        integrator_ida_pcreate_jacobian
252        , integrator_ida_psetup_jacobian
253        , integrator_ida_psolve_jacobian
254    };
255    
256    /*------
257    Jacobi preconditioner -- experimental    Jacobi preconditioner -- experimental
258  */  */
259    
# Line 1565  int integrator_ida_rootfn(realtype tt, N Line 1597  int integrator_ida_rootfn(realtype tt, N
1597                      gout[i] = -2.0;                      gout[i] = -2.0;
1598                  }                  }
1599                  break;                  break;
1600                case e_bnd_undefined:
1601                    ERROR_REPORTER_HERE(ASC_PROG_ERR,"Invalid boundary type e_bnd_undefined");
1602                    return 1;
1603          }          }
1604      }      }
1605    
1606      return 0; /* no way to detect errors in bndman_*_eval at this stage */      return 0; /* no way to detect errors in bndman_*_eval at this stage */
1607  }  }
1608    
1609    
1610    /*----------------------------------------------
1611      FULL JACOBIAN PRECONDITIONER -- EXPERIMENTAL.
1612    */
1613    
1614    void integrator_ida_pcreate_jacobian(IntegratorSystem *sys){
1615        IntegratorIdaData *enginedata =sys->enginedata;
1616        IntegratorIdaPrecDataJacobian *precdata;
1617        precdata = ASC_NEW(IntegratorIdaPrecDataJacobian);
1618        mtx_matrix_t P;
1619        asc_assert(sys->n_y);
1620        precdata->L = linsolqr_create_default();
1621    
1622        /* allocate matrix to be used by linsolqr */
1623        P = mtx_create();
1624        mtx_set_order(P, sys->n_y);
1625        linsolqr_set_matrix(precdata->L, P);
1626    
1627        enginedata->pfree = &integrator_ida_pfree_jacobian;
1628        enginedata->precdata = precdata;
1629        CONSOLE_DEBUG("Allocated memory for Full Jacobian preconditioner");
1630    }
1631    
1632    void integrator_ida_pfree_jacobian(IntegratorIdaData *enginedata){
1633        mtx_matrix_t P;
1634        IntegratorIdaPrecDataJacobian *precdata;
1635    
1636        if(enginedata->precdata){
1637            precdata = (IntegratorIdaPrecDataJacobian *)enginedata->precdata;
1638            P = linsolqr_get_matrix(precdata->L);
1639            mtx_destroy(P);
1640            linsolqr_destroy(precdata->L);
1641            ASC_FREE(precdata);
1642            enginedata->precdata = NULL;
1643    
1644            CONSOLE_DEBUG("Freed memory for Full Jacobian preconditioner");
1645        }
1646        enginedata->pfree = NULL;
1647    }
1648    
1649    /**
1650        EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1651    
1652        'setup' function.
1653    */
1654    int integrator_ida_psetup_jacobian(realtype tt,
1655             N_Vector yy, N_Vector yp, N_Vector rr,
1656             realtype c_j, void *p_data,
1657             N_Vector tmp1, N_Vector tmp2,
1658             N_Vector tmp3
1659    ){
1660        int i, j, res;
1661        IntegratorSystem *sys;
1662        IntegratorIdaData *enginedata;
1663        IntegratorIdaPrecDataJacobian *precdata;
1664        linsolqr_system_t L;
1665        mtx_matrix_t P;
1666    
1667        L = precdata->L;
1668        P = linsolqr_get_matrix(L);
1669        mtx_clear(P);
1670    
1671        struct rel_relation **relptr;
1672    
1673        sys = (IntegratorSystem *)p_data;
1674        enginedata = sys->enginedata;
1675        precdata = (IntegratorIdaPrecDataJacobian *)(enginedata->precdata);
1676        double *derivatives;
1677        struct var_variable **variables;
1678        int count, status;
1679        char *relname;
1680        mtx_coord_t C;
1681    
1682        CONSOLE_DEBUG("Setting up Jacobian preconditioner");
1683    
1684        variables = ASC_NEW_ARRAY(struct var_variable*, NV_LENGTH_S(yy) * 2);
1685        derivatives = ASC_NEW_ARRAY(double, NV_LENGTH_S(yy) * 2);
1686    
1687        /**
1688            @TODO FIXME here we are using the very inefficient and contorted approach
1689            of calculating the whole jacobian, then extracting just the diagonal elements.
1690        */
1691    
1692        for(i=0, relptr = enginedata->rellist;
1693                i< enginedata->nrels && relptr != NULL;
1694                ++i, ++relptr
1695        ){
1696            /* get derivatives for this particular relation */
1697            status = relman_diff3(*relptr, &enginedata->vfilter, derivatives, variables, &count, enginedata->safeeval);
1698            if(status){
1699                relname = rel_make_name(sys->system, *relptr);
1700                CONSOLE_DEBUG("ERROR calculating preconditioner derivatives for relation '%s'",relname);
1701                ASC_FREE(relname);
1702                break;
1703            }
1704            /* CONSOLE_DEBUG("Got %d derivatives from relation %d",count,i); */
1705            /* find the diagonal elements */
1706            for(j=0; j<count; ++j){
1707                if(var_deriv(variables[j])){
1708                    mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), c_j * derivatives[j]);
1709                }else{
1710                    mtx_fill_value(P, mtx_coord(&C, i, var_sindex(variables[j])), derivatives[j]);
1711                }
1712            }
1713        }
1714    
1715        mtx_assemble(P);
1716    
1717        if(status){
1718            CONSOLE_DEBUG("Error found when evaluating derivatives");
1719            res = 1; goto finish; /* recoverable */
1720        }
1721    
1722        integrator_ida_write_incidence(sys);
1723    
1724        res = 0;
1725    finish:
1726        ASC_FREE(variables);
1727        ASC_FREE(derivatives);
1728        return res;
1729    };
1730    
1731    /**
1732        EXPERIMENTAL. Full Jacobian preconditioner for use with IDA Krylov solvers
1733    
1734        'solve' function.
1735    */
1736    int integrator_ida_psolve_jacobian(realtype tt,
1737             N_Vector yy, N_Vector yp, N_Vector rr,
1738             N_Vector rvec, N_Vector zvec,
1739             realtype c_j, realtype delta, void *p_data,
1740             N_Vector tmp
1741    ){
1742        IntegratorSystem *sys;
1743        IntegratorIdaData *data;
1744        IntegratorIdaPrecDataJacobian *precdata;
1745        sys = (IntegratorSystem *)p_data;
1746        data = sys->enginedata;
1747        precdata = (IntegratorIdaPrecDataJacobian *)(data->precdata);
1748        linsolqr_system_t L = precdata->L;
1749    
1750        linsolqr_add_rhs(L,NV_DATA_S(rvec),FALSE);
1751    
1752        mtx_region_t R;
1753        R.row.low = R.col.low = 0;
1754        R.row.high = R.col.high = mtx_order(linsolqr_get_matrix(L)) - 1;
1755        linsolqr_set_region(L,R);
1756    
1757        linsolqr_prep(L,linsolqr_fmethod_to_fclass(linsolqr_fmethod(L)));
1758        linsolqr_reorder(L, &R, linsolqr_rmethod(L));      
1759    
1760        /// @TODO more here
1761    
1762        linsolqr_remove_rhs(L,NV_DATA_S(rvec));
1763    
1764        CONSOLE_DEBUG("Solving Jacobian preconditioner (c_j = %f)",c_j);
1765        return 0;
1766    };
1767    
1768    
1769  /*----------------------------------------------  /*----------------------------------------------
1770    JACOBI PRECONDITIONER -- EXPERIMENTAL.    JACOBI PRECONDITIONER -- EXPERIMENTAL.
1771  */  */
Legend:
Removed from v.1394  
changed lines
  Added in v.1395
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22