Changeset 9597


Ignore:
Timestamp:
Jan 27, 2017 4:52:26 PM (3 years ago)
Author:
Christian Andersson
Message:

Minor cleanup in runtime. Related to ticket:4815

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • PyFMI/trunk/src/pyfmi/simulation/assimulo_interface.py

    r9097 r9597  
    13711371
    13721372                f.write("Solver: %s \n"%solver.__class__.__name__)
    1373                 #f.write("State variables: "+names+ "\n")
     1373                f.write("State variables: "+names+ "\n")
    13741374
    13751375                str_y = ""
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c

    r9581 r9597  
    138138    }
    139139    return(sum);
     140}
     141
     142/* compare analytical/external and finite differences Jacobians */
     143static int jmi_kinsol_check_jacobian(jmi_block_solver_t* block, jmi_real_t* jac_finite_difference, jmi_real_t* jac_analytical) {
     144    int N = block->n, i, j;
     145    jmi_kinsol_solver_t* solver = block->solver;
     146   
     147    if (   block->dF
     148        ||(block->Jacobian && (solver->has_compression_setup_flag
     149        || block->options->jacobian_calculation_mode == jmi_calculate_externally_jacobian_calculation_mode)))
     150    {
     151        for (i = 0; i < N; i++) {
     152            for (j = 0; j < N; j++) {
     153                realtype fd_val = jac_finite_difference[i * N + j];
     154                realtype a_val = jac_analytical[i * N + j];
     155                realtype rel_error = RAbs(a_val - fd_val) / (RAbs(fd_val) + 1);
     156                if (rel_error >= block->options->block_jacobian_check_tol) {
     157                    if(block->options->jacobian_calculation_mode == jmi_calculate_externally_jacobian_calculation_mode) {
     158                        jmi_log_node(block->log, logError, "JacobianCheck",
     159                            "<j: %d, i: %d, external: %e, finiteDifference: %e, relativeError: %e>",
     160                            j, i, a_val, fd_val, rel_error);
     161                    }
     162                    else {
     163                        jmi_log_node(block->log, logError, "JacobianCheck",
     164                            "<j: %d, i: %d, analytic: %e, finiteDifference: %e, relativeError: %e>",
     165                            j, i, a_val, fd_val, rel_error);
     166                    }
     167                }
     168            }
     169        }
     170    } else {
     171        jmi_log_node(block->log, logError, "JacobianCheck",
     172            "No block jacobian specified, unable to do jacobian check");
     173    }
     174    return 0;
    140175}
    141176
     
    656691        }       
    657692    }
    658 
     693   
     694    /* Verify Jacobian */
    659695    if (block->options->block_jacobian_check) {
    660         /* compare analytical/external and finite differences Jacobians */
    661         if (block->dF || (block->Jacobian && (solver->has_compression_setup_flag
    662             || (block->options->jacobian_calculation_mode == jmi_calculate_externally_jacobian_calculation_mode)))) {
    663                 for (i = 0; i < N; i++) {
    664                     for (j = 0; j < N; j++) {
    665                         realtype fd_val = jac_fd[i * N + j];
    666                         realtype a_val = J->data[i * N + j];
    667                         realtype rel_error = RAbs(a_val - fd_val) / (RAbs(fd_val) + 1);
    668                         if (rel_error >= block->options->block_jacobian_check_tol) {
    669                             if(block->options->jacobian_calculation_mode == jmi_calculate_externally_jacobian_calculation_mode) {
    670                                 jmi_log_node(block->log, logError, "JacobianCheck",
    671                                     "<j: %d, i: %d, external: %e, finiteDifference: %e, relativeError: %e>",
    672                                     j, i, a_val, fd_val, rel_error);
    673                             }
    674                             else {
    675                                 jmi_log_node(block->log, logError, "JacobianCheck",
    676                                     "<j: %d, i: %d, analytic: %e, finiteDifference: %e, relativeError: %e>",
    677                                     j, i, a_val, fd_val, rel_error);
    678                             }
    679                         }
    680                     }
    681                 }
    682         } else {
    683             jmi_log_node(block->log, logError, "JacobianCheck",
    684                 "No block jacobian specified, unable to do jacobian check");
    685         }
     696        jmi_kinsol_check_jacobian(block, jac_fd, J->data);
    686697        free(jac_fd);
    687698    }
     
    12261237    KINSetMaxNewtonStep(solver->kin_mem, solver->max_nw_step);
    12271238   
    1228     if(block->options->iteration_variable_scaling_mode)
    1229     {
    1230         /*
    1231             Set variable scaling based on nominal values.         
    1232         */
     1239    /* Set variable scaling based on nominal values. */
     1240    if(block->options->iteration_variable_scaling_mode) {
    12331241        int i;
    12341242        for(i=0;i< block->n;++i){
     
    17181726    double J_recip_cond = 1.0;
    17191727    int info;
    1720 
     1728   
     1729    /* Compute infinity norm of J to be used with dgecon */
     1730    J_norm = dlange_(&norm, &N, &N, block->J->data, &N, solver->lapack_work);
     1731   
    17211732    /* Copy Jacobian to factorization matrix */
    17221733    DenseCopy(block->J, solver->J_LU);
     
    17271738        return 1e100;
    17281739    }
    1729 
    1730     /* Compute infinity norm of J to be used with dgecon */
    1731     J_norm = dlange_(&norm, &N, &N, block->J->data, &N, solver->lapack_work);
    1732 
    17331740    /* Compute reciprocal condition number */
    17341741    dgecon_(&norm, &N, solver->J_LU->data, &N, &J_norm, &J_recip_cond, solver->lapack_work, solver->lapack_iwork,&info);
     
    24322439    solver->J_LU = NewDenseMat(n ,n);
    24332440    solver->J_sing = NewDenseMat(n, n);
    2434     solver->J_SVD_U = NewDenseMat(n, n);
    2435     solver->J_SVD_VT = NewDenseMat(n, n);
    24362441    solver->J_Dependency = NewDenseMat(n,n);
    24372442    solver->J_is_singular_flag = 0;
     
    24722477   
    24732478    /*Attach linear solver*/
    2474     /*Dense Kinsol solver*/
    2475     /*flag = KINDense(solver->kin_mem, block->n);
    2476       jmi_kinsol_error_handling(flag);*/
    2477      
    2478      
    2479    /*Dense Kinsol using regularization*/
    2480     /*    flag = KINPinv(solver->kin_mem, block->n);
    2481     jmi_kinsol_error_handling(jmi, flag);
    2482     KINDlsSetDenseJacFn(solver->kin_mem, (KINDlsDenseJacFn)kin_dF);
    2483 
    2484     */
    24852479    kin_mem->kin_lsetup = jmi_kin_lsetup;
    24862480    kin_mem->kin_lsolve = jmi_kin_lsolve;
    24872481    kin_mem->kin_setupNonNull = TRUE;
    24882482    kin_mem->kin_inexact_ls = FALSE;
    2489    
    24902483    /*End linear solver*/
    24912484   
     
    25572550    DestroyMat(solver->J_LU);
    25582551    DestroyMat(solver->J_sing);
    2559     DestroyMat(solver->J_SVD_U);
    2560     DestroyMat(solver->J_SVD_VT);
    25612552    DestroyMat(solver->J_Dependency);
    25622553    free(solver->cScale);
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

    r9581 r9597  
    9898    DlsMat J_LU;                    /**< \brief Jacobian matrix/it's LU decomposition */
    9999    DlsMat J_sing;                  /**< \brief Jacobian matrix/it's right singular vectors */
    100     DlsMat J_SVD_U;                 /**< \brief The left singular vectors */
    101     DlsMat J_SVD_VT;                /**< \brief The right singular vectors */
    102100    DlsMat J_Dependency;            /**< \brief Dependency matrix with value 1 at (i,j) if iv j depends on residual i, 0 otherwise */
    103101
  • trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.c

    r9571 r9597  
    378378        }
    379379        info = block->F(block->problem_data,block->x, solver->rhs, JMI_BLOCK_EVALUATE);
    380         /* info = block->F(block->problem_data,block->last_accepted_x, solver->rhs, JMI_BLOCK_EVALUATE); */
    381380    } else {
    382381        /* Ignore bounds when calculating RHS with zero vector*/
Note: See TracChangeset for help on using the changeset viewer.