From e7d386dfb5d12bf8923cc76196e78a0a3048dfcc Mon Sep 17 00:00:00 2001 From: Willi Braun Date: Wed, 5 Dec 2018 12:31:48 +0100 Subject: [PATCH] [Backend] symbolic jacobian remove defines in generates code - prepare for thread-safe linear systems in symbolic jacobian - mark which linear systems are part of the jacobian - expand most jacobian related functions in the cRuntime Belonging to [master]: - OpenModelica/OMCompiler#2745 - OpenModelica/OpenModelica-testsuite#1084 --- Compiler/BackEnd/BackendDAEUtil.mo | 10 + Compiler/BackEnd/Differentiate.mo | 37 +- Compiler/BackEnd/HpcOmScheduler.mo | 5 +- Compiler/BackEnd/SymbolicJacobian.mo | 63 ++- Compiler/FrontEnd/ComponentReference.mo | 31 ++ Compiler/FrontEnd/DAE.mo | 1 + Compiler/SimCode/ReduceDAE.mo | 7 +- Compiler/SimCode/SimCode.mo | 4 +- Compiler/SimCode/SimCodeFunction.mo | 5 +- Compiler/SimCode/SimCodeUtil.mo | 148 ++++-- Compiler/Stubs/SimCodeUtil.mo | 8 + Compiler/Template/CodegenC.tpl | 502 ++++++++++-------- Compiler/Template/CodegenCFunctions.tpl | 47 +- Compiler/Template/SimCodeTV.mo | 56 ++ Compiler/Util/Error.mo | 3 + .../c/linearization/linearize.cpp | 117 ++-- SimulationRuntime/c/openmodelica_func.h | 20 +- .../DataManagement/DerStructure.c | 8 +- .../c/optimization/DataManagement/MoveData.c | 36 +- SimulationRuntime/c/simulation/solver/dassl.c | 68 +-- .../c/simulation/solver/ida_solver.c | 7 +- .../c/simulation/solver/kinsolSolver.c | 2 +- .../c/simulation/solver/linearSolverKlu.c | 39 +- .../c/simulation/solver/linearSolverKlu.h | 2 +- .../c/simulation/solver/linearSolverLapack.c | 35 +- .../c/simulation/solver/linearSolverLapack.h | 2 +- .../c/simulation/solver/linearSolverLis.c | 42 +- .../c/simulation/solver/linearSolverLis.h | 2 +- .../solver/linearSolverTotalPivot.c | 44 +- .../solver/linearSolverTotalPivot.h | 2 +- .../c/simulation/solver/linearSolverUmfpack.c | 51 +- .../c/simulation/solver/linearSolverUmfpack.h | 2 +- .../c/simulation/solver/linearSystem.c | 29 +- .../c/simulation/solver/linearSystem.h | 2 +- .../solver/nonlinearSolverHomotopy.c | 29 +- .../simulation/solver/nonlinearSolverHybrd.c | 29 +- .../simulation/solver/nonlinearSolverNewton.c | 29 +- .../c/simulation/solver/nonlinearSystem.c | 3 +- .../solver/perform_qss_simulation.c | 4 +- .../c/simulation/solver/stateset.c | 60 ++- SimulationRuntime/c/simulation_data.h | 14 +- .../fmi/export/fmi2/fmu2_model_interface.c | 5 +- 42 files changed, 962 insertions(+), 648 deletions(-) diff --git a/Compiler/BackEnd/BackendDAEUtil.mo b/Compiler/BackEnd/BackendDAEUtil.mo index d5bc5e908b9..c98b1fe562d 100644 --- a/Compiler/BackEnd/BackendDAEUtil.mo +++ b/Compiler/BackEnd/BackendDAEUtil.mo @@ -142,6 +142,16 @@ algorithm end match; end isSimulationDAE; +public function isJacobianDAE + input BackendDAE.Shared inShared; + output Boolean res; +algorithm + res := match(inShared) + case (BackendDAE.SHARED(backendDAEType=BackendDAE.JACOBIAN())) then true; + else false; + end match; +end isJacobianDAE; + /************************************************* * checkBackendDAE and stuff ************************************************/ diff --git a/Compiler/BackEnd/Differentiate.mo b/Compiler/BackEnd/Differentiate.mo index c1cefd9bf26..557c9a2403e 100644 --- a/Compiler/BackEnd/Differentiate.mo +++ b/Compiler/BackEnd/Differentiate.mo @@ -1200,7 +1200,7 @@ algorithm //Take care! state means => der(state) false = BackendVariable.isStateVar(var); - cr = createDifferentiatedCrefName(cr, inDiffwrtCref, matrixName); + cr = ComponentReference.createDifferentiatedCrefName(cr, inDiffwrtCref, matrixName); res = DAE.CREF(cr, tp); then (res, inFunctionTree); @@ -1212,7 +1212,7 @@ algorithm //Take care! state means => der(state) false = BackendVariable.isStateVar(var); - cr = createDifferentiatedCrefName(cr, inDiffwrtCref, matrixName); + cr = ComponentReference.createDifferentiatedCrefName(cr, inDiffwrtCref, matrixName); res = DAE.CREF(cr, tp); then (res, inFunctionTree); @@ -1257,37 +1257,6 @@ algorithm outCref := ComponentReference.crefSetLastType(outCref, ComponentReference.crefLastType(inCref)); end createDiffedCrefName; -public function createDifferentiatedCrefName - input DAE.ComponentRef inCref; - input DAE.ComponentRef inX; - input String inMatrixName; - output DAE.ComponentRef outCref; -protected - list subs; - constant Boolean debug = false; -algorithm - if debug then print("inCref: " + ComponentReference.printComponentRefStr(inCref) +"\n"); end if; - - // move subs and and type to lastCref, to move type replace by last type - // and move last cref type to the last cref. - subs := ComponentReference.crefLastSubs(inCref); - outCref := ComponentReference.crefStripLastSubs(inCref); - outCref := ComponentReference.replaceSubsWithString(outCref); - if debug then print("after full type " + Types.printTypeStr(ComponentReference.crefTypeConsiderSubs(inCref)) + "\n"); end if; - outCref := ComponentReference.crefSetLastType(outCref, DAE.T_UNKNOWN_DEFAULT); - if debug then print("after strip: " + ComponentReference.printComponentRefListStr(ComponentReference.expandCref(outCref, true)) + "\n"); end if; - - // join crefs - outCref := ComponentReference.joinCrefs(outCref, ComponentReference.makeCrefIdent(BackendDAE.partialDerivativeNamePrefix + inMatrixName, DAE.T_UNKNOWN_DEFAULT, {})); - outCref := ComponentReference.joinCrefs(outCref, inX); - if debug then print("after join: " + ComponentReference.printComponentRefListStr(ComponentReference.expandCref(outCref, true)) + "\n"); end if; - - // fix subs and type of the last cref - outCref := ComponentReference.crefSetLastSubs(outCref, subs); - outCref := ComponentReference.crefSetLastType(outCref, ComponentReference.crefLastType(inCref)); - if debug then print("outCref: " + ComponentReference.printComponentRefStr(outCref) +"\n"); end if; -end createDifferentiatedCrefName; - public function createSeedCrefName input DAE.ComponentRef inCref; input String inMatrixName; @@ -1382,7 +1351,7 @@ algorithm cr = Expression.expCref(e); tp = Expression.typeof(e); cr = ComponentReference.crefPrefixDer(cr); - cr = createDifferentiatedCrefName(cr, inDiffwrtCref, matrixName); + cr = ComponentReference.createDifferentiatedCrefName(cr, inDiffwrtCref, matrixName); res = Expression.makeCrefExp(cr, tp); b = ComponentReference.crefEqual(DAE.CREF_IDENT("$",DAE.T_REAL_DEFAULT,{}), inDiffwrtCref); diff --git a/Compiler/BackEnd/HpcOmScheduler.mo b/Compiler/BackEnd/HpcOmScheduler.mo index 515b9753ac0..d0b01d4b2aa 100644 --- a/Compiler/BackEnd/HpcOmScheduler.mo +++ b/Compiler/BackEnd/HpcOmScheduler.mo @@ -3550,10 +3550,11 @@ algorithm list> colCols; array ass; Integer newIdx; - case(SOME(SimCode.JAC_MATRIX(jacCols,vars,name,sparsity,sparsityT,colCols,maxCol,jacIdx,partIdx)),(newIdx,ass)) + Option crefToSimVarHTJacobian; + case(SOME(SimCode.JAC_MATRIX(jacCols,vars,name,sparsity,sparsityT,colCols,maxCol,jacIdx,partIdx,crefToSimVarHTJacobian)),(newIdx,ass)) equation (jacCols,(newIdx,ass)) = List.mapFold(jacCols,TDS_replaceSimEqSysIdxInJacobianColumnWithUpdate,(newIdx,ass)); - then (SOME(SimCode.JAC_MATRIX(jacCols,vars,name,sparsity,sparsityT,colCols,maxCol,jacIdx,partIdx)),(newIdx,ass)); + then (SOME(SimCode.JAC_MATRIX(jacCols,vars,name,sparsity,sparsityT,colCols,maxCol,jacIdx,partIdx,crefToSimVarHTJacobian)),(newIdx,ass)); else (jacIn,tplIn); end matchcontinue; end TDS_replaceSimEqSysIdxInJacobianMatrixWithUpdate; diff --git a/Compiler/BackEnd/SymbolicJacobian.mo b/Compiler/BackEnd/SymbolicJacobian.mo index 59180a40ca0..d4a0ac67abb 100644 --- a/Compiler/BackEnd/SymbolicJacobian.mo +++ b/Compiler/BackEnd/SymbolicJacobian.mo @@ -1972,6 +1972,7 @@ algorithm outFunctionTree := shared.functionTree; if not onlySparsePattern then (symbolicJacobian, outFunctionTree) := createJacobian(inBackendDAE,inDiffVars, inStateVars, inInputVars, inParameterVars, inDifferentiatedVars, inVars, inName); + true := checkForNonLinearStrongComponents(symbolicJacobian); outJacobian := SOME(symbolicJacobian); else outJacobian := NONE(); @@ -2294,14 +2295,14 @@ algorithm try (_, _) := BackendVariable.getVarSingle(currVar, inAllVars); currVar := ComponentReference.crefPrefixDer(currVar); - derivedCref := Differentiate.createDifferentiatedCrefName(currVar, cref, inMatrixName); + derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName); r1 := BackendVariable.copyVarNewName(derivedCref, v); r1 := BackendVariable.setVarKind(r1, BackendDAE.STATE_DER()); r1.unreplaceable := true; index := index + 1; else currVar := ComponentReference.crefPrefixDer(currVar); - derivedCref := Differentiate.createDifferentiatedCrefName(currVar, cref, inMatrixName); + derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName); r1 := BackendVariable.copyVarNewName(derivedCref, v); r1 := BackendVariable.setVarKind(r1, BackendDAE.STATE_DER()); end try; @@ -2311,13 +2312,13 @@ algorithm case((v as BackendDAE.VAR(varName=currVar))::restVar, cref, _, index, _, _) algorithm try (_, _) := BackendVariable.getVarSingle(currVar, inAllVars); - derivedCref := Differentiate.createDifferentiatedCrefName(currVar, cref, inMatrixName); + derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName); r1 := BackendVariable.copyVarNewName(derivedCref, v); r1 := BackendVariable.setVarKind(r1, BackendDAE.VARIABLE()); r1.unreplaceable := true; index := index + 1; else - derivedCref := Differentiate.createDifferentiatedCrefName(currVar, cref, inMatrixName); + derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName); r1 := BackendVariable.copyVarNewName(derivedCref, v); r1 := BackendVariable.setVarKind(r1, BackendDAE.VARIABLE()); end try; @@ -3700,5 +3701,59 @@ algorithm end matchcontinue; end rhsConstant2; +// ============================================================================= +// Function detects non-linear strong component in symbolic jacobians +// - non-linear components should never appear in symbolic jacobian and +// indicate an singular or wrong system +// - this modules stops compiling and outputs an error, otherwise we +// would get error at runtime compiling +// ============================================================================= + +function checkForNonLinearStrongComponents +"Checks for non-linear algebraic strong compontents and break if some found." + input BackendDAE.SymbolicJacobian symbolicJacobian; + output Boolean result; +protected + BackendDAE.BackendDAE jacBDAE; + String name; +algorithm + (jacBDAE, name, _, _, _) := symbolicJacobian; + try + _ := BackendDAEUtil.mapEqSystem(jacBDAE, checkForNonLinearStrongComponents_work); + result := true; + else + Error.addMessage(Error.INVALID_NONLINEAR_JACOBIAN_COMPONENT, {name}); + result := false; + end try; +end checkForNonLinearStrongComponents; + +function checkForNonLinearStrongComponents_work + input output BackendDAE.EqSystem syst; + input output BackendDAE.Shared shared; +protected + BackendDAE.StrongComponents comps; +algorithm + try + BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(comps=comps)) := syst; + for comp in comps loop + () := match (comp) + local + BackendDAE.JacobianType jacTp; + case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_NONLINEAR()) + then fail(); + case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_NO_ANALYTIC()) + then fail(); + case BackendDAE.EQUATIONSYSTEM(jacType=BackendDAE.JAC_GENERIC()) + then fail(); + case BackendDAE.TORNSYSTEM(linear=false) + then fail(); + else (); + end match; + end for; + else + fail(); + end try; +end checkForNonLinearStrongComponents_work; + annotation(__OpenModelica_Interface="backend"); end SymbolicJacobian; diff --git a/Compiler/FrontEnd/ComponentReference.mo b/Compiler/FrontEnd/ComponentReference.mo index 7858995f722..9e45f266a51 100644 --- a/Compiler/FrontEnd/ComponentReference.mo +++ b/Compiler/FrontEnd/ComponentReference.mo @@ -3955,5 +3955,36 @@ algorithm end while; end getConsumedMemory; +public function createDifferentiatedCrefName + input DAE.ComponentRef inCref; + input DAE.ComponentRef inX; + input String inMatrixName; + output DAE.ComponentRef outCref; +protected + list subs; + constant Boolean debug = false; +algorithm + if debug then print("inCref: " + ComponentReference.printComponentRefStr(inCref) +"\n"); end if; + + // move subs and and type to lastCref, to move type replace by last type + // and move last cref type to the last cref. + subs := ComponentReference.crefLastSubs(inCref); + outCref := ComponentReference.crefStripLastSubs(inCref); + outCref := ComponentReference.replaceSubsWithString(outCref); + if debug then print("after full type " + Types.printTypeStr(ComponentReference.crefTypeConsiderSubs(inCref)) + "\n"); end if; + outCref := ComponentReference.crefSetLastType(outCref, DAE.T_UNKNOWN_DEFAULT); + if debug then print("after strip: " + ComponentReference.printComponentRefListStr(ComponentReference.expandCref(outCref, true)) + "\n"); end if; + + // join crefs + outCref := ComponentReference.joinCrefs(outCref, ComponentReference.makeCrefIdent(DAE.partialDerivativeNamePrefix + inMatrixName, DAE.T_UNKNOWN_DEFAULT, {})); + outCref := ComponentReference.joinCrefs(outCref, inX); + if debug then print("after join: " + ComponentReference.printComponentRefListStr(ComponentReference.expandCref(outCref, true)) + "\n"); end if; + + // fix subs and type of the last cref + outCref := ComponentReference.crefSetLastSubs(outCref, subs); + outCref := ComponentReference.crefSetLastType(outCref, ComponentReference.crefLastType(inCref)); + if debug then print("outCref: " + ComponentReference.printComponentRefStr(outCref) +"\n"); end if; +end createDifferentiatedCrefName; + annotation(__OpenModelica_Interface="frontend"); end ComponentReference; diff --git a/Compiler/FrontEnd/DAE.mo b/Compiler/FrontEnd/DAE.mo index d022d333df5..dc060b3ca72 100644 --- a/Compiler/FrontEnd/DAE.mo +++ b/Compiler/FrontEnd/DAE.mo @@ -61,6 +61,7 @@ public type StartValue = Option; public constant String UNIQUEIO = "$unique$outer$"; public constant String derivativeNamePrefix = "$DER"; +public constant String partialDerivativeNamePrefix = "$pDER"; public constant String preNamePrefix = "$PRE"; public constant String previousNamePrefix = "$CLKPRE"; public constant String startNamePrefix = "$START"; diff --git a/Compiler/SimCode/ReduceDAE.mo b/Compiler/SimCode/ReduceDAE.mo index 83e97f54322..1d3dae9607b 100644 --- a/Compiler/SimCode/ReduceDAE.mo +++ b/Compiler/SimCode/ReduceDAE.mo @@ -342,6 +342,7 @@ protected function addLabelToEquations list residual; Option jacobianMatrix; BackendDAE.EquationAttributes eqAttr; + Boolean partOfJac; // nothing case ({},vars,idx,_,_) then ({},vars,idx,{}); @@ -395,8 +396,9 @@ protected function addLabelToEquations labels3=listAppend(labels,labels2); then (SimCode.SES_ALGORITHM(i,statements2, eqAttr) :: es_1,vars_2,idx3,labels3); + // linear systems - case (((eq as SimCode.SES_LINEAR (SimCode.LINEARSYSTEM(i,partOfLinear,tornSystem,varsLin,b,A,residual,jacobianMatrix,sourcelist,idxLS,nUnknownsLS),NONE(), eqAttr)) :: es),vars,idx,_,_) + case (((eq as SimCode.SES_LINEAR (SimCode.LINEARSYSTEM(i,partOfLinear,tornSystem,varsLin,b,A,residual,jacobianMatrix,sourcelist,idxLS,nUnknownsLS,partOfJac),NONE(), eqAttr)) :: es),vars,idx,_,_) equation if(Flags.isSet(Flags.REDUCE_DAE)) then @@ -409,7 +411,8 @@ protected function addLabelToEquations labels3=listAppend(labels,labels2); then - (SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(i,partOfLinear,tornSystem,varsLin,b,A2,residual,jacobianMatrix,sourcelist,idxLS,nUnknownsLS),NONE(), eqAttr) :: es_1,vars_2,idx3,labels3); + (SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(i,partOfLinear,tornSystem,varsLin,b,A2,residual,jacobianMatrix,sourcelist,idxLS,nUnknownsLS, partOfJac),NONE(), eqAttr) :: es_1,vars_2,idx3,labels3); + // non-linear systems case (((eq as SimCode.SES_NONLINEAR(SimCode.NONLINEARSYSTEM(index=i,eqs=nl,crefs=crefs,indexNonLinearSystem=idxNLS,nUnknowns=nUnknownsNLS,jacobianMatrix=jacobianMatrix),NONE(), eqAttr)) :: es),vars,idx,_,_) equation diff --git a/Compiler/SimCode/SimCode.mo b/Compiler/SimCode/SimCode.mo index f5bb947e207..8a867bdb9b0 100644 --- a/Compiler/SimCode/SimCode.mo +++ b/Compiler/SimCode/SimCode.mo @@ -88,10 +88,11 @@ uniontype JacobianMatrix Integer maxColorCols; Integer jacobianIndex; Integer partitionIndex; + Option crefsHT; // all jacobian variables end JAC_MATRIX; end JacobianMatrix; -constant JacobianMatrix emptyJacobian = JAC_MATRIX({}, {}, "", {}, {}, {}, 0, -1, 0); +constant JacobianMatrix emptyJacobian = JAC_MATRIX({}, {}, "", {}, {}, {}, 0, -1, 0, NONE()); constant PartitionData emptyPartitionData = PARTITIONDATA(-1,{},{},{}); @@ -411,6 +412,7 @@ uniontype LinearSystem list sources; Integer indexLinearSystem; Integer nUnknowns "Number of variables that are solved in this system. Needed because 'crefs' only contains the iteration variables."; + Boolean partOfJac "if TRUE then this system is part of a jacobian matrix"; end LINEARSYSTEM; end LinearSystem; diff --git a/Compiler/SimCode/SimCodeFunction.mo b/Compiler/SimCode/SimCodeFunction.mo index e41f1a6ad5a..6683345c840 100644 --- a/Compiler/SimCode/SimCodeFunction.mo +++ b/Compiler/SimCode/SimCodeFunction.mo @@ -37,6 +37,7 @@ public import Absyn; import DAE; import HashTableStringToPath; +import HashTableCrefSimVar; import Tpl; uniontype FunctionCode @@ -206,6 +207,7 @@ uniontype Context end ALGLOOP_CONTEXT; record JACOBIAN_CONTEXT + Option jacHT; end JACOBIAN_CONTEXT; record OTHER_CONTEXT @@ -225,12 +227,13 @@ uniontype Context record DAE_MODE_CONTEXT end DAE_MODE_CONTEXT; + end Context; public constant Context contextSimulationNonDiscrete = SIMULATION_CONTEXT(false); public constant Context contextSimulationDiscrete = SIMULATION_CONTEXT(true); public constant Context contextFunction = FUNCTION_CONTEXT(); -public constant Context contextJacobian = JACOBIAN_CONTEXT(); +public constant Context contextJacobian = JACOBIAN_CONTEXT(NONE()); public constant Context contextAlgloopJacobian = ALGLOOP_CONTEXT(false,true); public constant Context contextAlgloopInitialisation = ALGLOOP_CONTEXT(true,false); public constant Context contextAlgloop = ALGLOOP_CONTEXT(false,false); diff --git a/Compiler/SimCode/SimCodeUtil.mo b/Compiler/SimCode/SimCodeUtil.mo index 9adfdaf3c75..0a490c796ea 100644 --- a/Compiler/SimCode/SimCodeUtil.mo +++ b/Compiler/SimCode/SimCodeUtil.mo @@ -3147,6 +3147,7 @@ algorithm Boolean mixedSystem; BackendDAE.TearingSet strictTearingSet; Option casualTearingSet; + Boolean partOfJac; // EQUATIONSYSTEM: continuous system of equations case (BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns), @@ -3170,7 +3171,8 @@ algorithm var_lst_1 = List.map(var_lst, BackendVariable.transformXToXd); vars_1 = BackendVariable.listVar1(var_lst_1); eqns_1 = BackendEquation.listEquation(eqn_lst); - (equations_, uniqueEqIndex, tempvars) = createOdeSystem2(false, vars_1, globalKnownVars, eqns_1, jacobian, jac_tp, funcs, vars, iuniqueEqIndex, ei, mixedSystem, itempvars); + partOfJac = BackendDAEUtil.isJacobianDAE(ishared); + (equations_, uniqueEqIndex, tempvars) = createOdeSystem2(false, vars_1, globalKnownVars, eqns_1, jacobian, jac_tp, funcs, vars, iuniqueEqIndex, ei, mixedSystem, partOfJac, itempvars); uniqueEqIndexMapping = uniqueEqIndex-1; //a system with this index is created that contains all the equations with the indeces from iuniqueEqIndex to uniqueEqIndex-2 //tmpEqSccMapping = appendSccIdxRange(iuniqueEqIndex, uniqueEqIndex - 1, isccIndex, ieqSccMapping); tmpEqSccMapping = appendSccIdxRange(uniqueEqIndexMapping, uniqueEqIndex - 1, isccIndex, ieqSccMapping); @@ -3209,6 +3211,7 @@ protected function createOdeSystem2 input Integer iuniqueEqIndex; input BackendDAE.ExtraInfo iei; input Boolean mixedSystem; + input Boolean partOfJac; input list itempvars; output list equations_; output Integer ouniqueEqIndex; @@ -3265,11 +3268,11 @@ algorithm (beqs, sources) = BackendDAEUtil.getEqnSysRhs(inEquationArray, inVars, SOME(inFuncs)); beqs = listReverse(beqs); simJac = List.map1(jac, jacToSimjac, inVars); - if (Config.simCodeTarget()=="Cpp") then - simJac = List.sort(simJac,simJacCSRToCSC); - end if; + if (Config.simCodeTarget()=="Cpp") then + simJac = List.sort(simJac,simJacCSRToCSC); + end if; - then ({SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(iuniqueEqIndex, mixedEvent, false, simVars, beqs, simJac, {}, NONE(), sources, 0, inVars.numberOfVars), NONE(), BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN)}, iuniqueEqIndex+1, itempvars); + then ({SimCode.SES_LINEAR(SimCode.LINEARSYSTEM(iuniqueEqIndex, mixedEvent, false, simVars, beqs, simJac, {}, NONE(), sources, 0, inVars.numberOfVars, partOfJac), NONE(), BackendDAE.EQ_ATTR_DEFAULT_UNKNOWN)}, iuniqueEqIndex+1, itempvars); // Time varying nonlinear jacobian. Non-linear system of equations. case (_, BackendDAE.JAC_GENERIC()) equation @@ -3387,6 +3390,8 @@ algorithm SimCode.NonlinearSystem nlSystem; Option alternativeTearingL; Option alternativeTearingNl; + BackendDAE.BackendDAEType backendDAEType; + Boolean partOfJac; // CASE: linear case(true, BackendDAE.EQSYSTEM(orderedVars=vars, orderedEqs=eqns), BackendDAE.SHARED(globalKnownVars=globalKnownVars)) equation @@ -3412,7 +3417,8 @@ algorithm simequations = listAppend(simequations, resEqs); (jacobianMatrix, uniqueEqIndex, tempvars) = createSymbolicSimulationJacobian(inJacobian, uniqueEqIndex, tempvars); - lSystem = SimCode.LINEARSYSTEM(uniqueEqIndex, false, true, simVars, {}, {}, simequations, jacobianMatrix, {}, 0, listLength(tvars)+nInnerVars+listLength(tempvars)-listLength(itempvars)); + partOfJac = BackendDAEUtil.isJacobianDAE(ishared); + lSystem = SimCode.LINEARSYSTEM(uniqueEqIndex, false, true, simVars, {}, {}, simequations, jacobianMatrix, {}, 0, listLength(tvars)+nInnerVars+listLength(tempvars)-listLength(itempvars), partOfJac); tempvars2 = tempvars; // Do if dynamic tearing is activated @@ -3433,7 +3439,7 @@ algorithm simequations = listAppend(simequations, resEqs); (jacobianMatrix, uniqueEqIndex, tempvars2) = createSymbolicSimulationJacobian(inJacobian, uniqueEqIndex, tempvars2); - alternativeTearingL = SOME(SimCode.LINEARSYSTEM(uniqueEqIndex, false, true, simVars, {}, {}, simequations, jacobianMatrix, {}, 0, listLength(tvars)+nInnerVars+listLength(tempvars2)-listLength(tempvars))); + alternativeTearingL = SOME(SimCode.LINEARSYSTEM(uniqueEqIndex, false, true, simVars, {}, {}, simequations, jacobianMatrix, {}, 0, listLength(tvars)+nInnerVars+listLength(tempvars2)-listLength(tempvars), partOfJac)); else alternativeTearingL = NONE(); @@ -3895,7 +3901,7 @@ algorithm DAE.FunctionTree funcs; - SimCode.HashTableCrefToSimVar crefSimVarHT; + HashTableCrefSimVar.HashTable crefToSimVarHTJacobian; case (BackendDAE.EMPTY_JACOBIAN(), _, _) then (NONE(), iuniqueEqIndex, itempvars); @@ -3934,7 +3940,7 @@ algorithm print("created sparse pattern for algebraic loop time: " + realString(clock()) + "\n"); end if; - then (SOME(SimCode.JAC_MATRIX({}, {}, "", sparseInts, sparseIntsT, coloring, maxColor, -1, 0)), iuniqueEqIndex, itempvars); + then (SOME(SimCode.JAC_MATRIX({}, {}, "", sparseInts, sparseIntsT, coloring, maxColor, -1, 0, NONE())), iuniqueEqIndex, itempvars); case (BackendDAE.GENERIC_JACOBIAN(SOME((BackendDAE.DAE(eqs={syst as BackendDAE.EQSYSTEM(matching=BackendDAE.MATCHING(comps=comps))}, shared=shared), name, @@ -4005,7 +4011,11 @@ algorithm print("analytical Jacobians -> transformed to SimCode for Matrix " + name + " time: " + realString(clock()) + "\n"); end if; - then (SOME(SimCode.JAC_MATRIX({SimCode.JAC_COLUMN(columnEquations, columnVars, nRows)}, seedVars, name, sparseInts, sparseIntsT, coloring, maxColor, -1, 0)), uniqueEqIndex, tempvars); + crefToSimVarHTJacobian = HashTableCrefSimVar.emptyHashTableSized(listLength(seedVars)+ listLength(columnVars)); + crefToSimVarHTJacobian = List.fold(seedVars, addSimVarToHashTable, crefToSimVarHTJacobian); + crefToSimVarHTJacobian = List.fold(columnVars, addSimVarToHashTable, crefToSimVarHTJacobian); + + then (SOME(SimCode.JAC_MATRIX({SimCode.JAC_COLUMN(columnEquations, columnVars, nRows)}, seedVars, name, sparseInts, sparseIntsT, coloring, maxColor, -1, 0, SOME(crefToSimVarHTJacobian))), uniqueEqIndex, tempvars); else equation @@ -4119,6 +4129,8 @@ algorithm Option optionBDAE; SimCode.JacobianMatrix tmpJac; + HashTableCrefSimVar.HashTable crefToSimVarHTJacobian; + case (_, _, _, {}) then (inJacobianMatrixes, iuniqueEqIndex); // if nothing is generated @@ -4194,7 +4206,7 @@ algorithm seedVars = List.map1(seedVars, setSimVarKind, BackendDAE.SEED_VAR()); seedVars = List.map1(seedVars, setSimVarMatrixName, SOME(name)); - tmpJac = SimCode.JAC_MATRIX({SimCode.JAC_COLUMN({},{},nRows)}, seedVars, name, sparseInts, sparseIntsT, coloring, maxColor, -1, 0); + tmpJac = SimCode.JAC_MATRIX({SimCode.JAC_COLUMN({},{},nRows)}, seedVars, name, sparseInts, sparseIntsT, coloring, maxColor, -1, 0, NONE()); linearModelMatrices = tmpJac::inJacobianMatrixes; (linearModelMatrices, uniqueEqIndex) = createSymbolicJacobianssSimCode(rest, inSimVarHT, iuniqueEqIndex, restnames, linearModelMatrices); @@ -4295,7 +4307,12 @@ algorithm seedVars = List.map1(seedVars, setSimVarKind, BackendDAE.SEED_VAR()); seedVars = List.map1(seedVars, setSimVarMatrixName, SOME(name)); - tmpJac = SimCode.JAC_MATRIX({SimCode.JAC_COLUMN(columnEquations, columnVars, nRows)}, seedVars, name, sparseInts, sparseIntsT, coloring, maxColor, -1, 0); + // create hash table for this jacobians + crefToSimVarHTJacobian = HashTableCrefSimVar.emptyHashTableSized(listLength(seedVars)+ listLength(columnVars)); + crefToSimVarHTJacobian = List.fold(seedVars, addSimVarToHashTable, crefToSimVarHTJacobian); + crefToSimVarHTJacobian = List.fold(columnVars, addSimVarToHashTable, crefToSimVarHTJacobian); + + tmpJac = SimCode.JAC_MATRIX({SimCode.JAC_COLUMN(columnEquations, columnVars, nRows)}, seedVars, name, sparseInts, sparseIntsT, coloring, maxColor, -1, 0, SOME(crefToSimVarHTJacobian)); linearModelMatrices = tmpJac::inJacobianMatrixes; (linearModelMatrices, uniqueEqIndex) = createSymbolicJacobianssSimCode(rest, inSimVarHT, uniqueEqIndex, restnames, linearModelMatrices); then @@ -4383,7 +4400,7 @@ algorithm case BackendDAE.STATE() then ComponentReference.crefPrefixDer(currVar); else currVar; end match; - derivedCref := Differentiate.createDifferentiatedCrefName(currVar, cref, inMatrixName); + derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName); v1 := BackendVariable.copyVarNewName(derivedCref, v); v1 := BackendVariable.setVarKind(v1, BackendDAE.JAC_VAR()); simVar := dlowvarToSimvar(v1, NONE(), inAllVars); @@ -4394,7 +4411,7 @@ algorithm case BackendDAE.STATE() then ComponentReference.crefPrefixDer(currVar); else currVar; end match; - derivedCref := Differentiate.createDifferentiatedCrefName(currVar, cref, inMatrixName); + derivedCref := ComponentReference.createDifferentiatedCrefName(currVar, cref, inMatrixName); v1 := BackendVariable.copyVarNewName(derivedCref, v); v1 := BackendVariable.setVarKind(v1, BackendDAE.JAC_DIFF_VAR()); simVar := dlowvarToSimvar(v1, NONE(), inAllVars); @@ -12666,50 +12683,71 @@ public function cref2simvar input DAE.ComponentRef inCref; input SimCode.SimCode simCode; output SimCodeVar.SimVar outSimVar; +protected + HashTableCrefSimVar.HashTable crefToSimVarHT; + DAE.ComponentRef badcref; algorithm - outSimVar := matchcontinue(inCref, simCode) - local - DAE.ComponentRef cref, badcref; - SimCodeVar.SimVar sv; - SimCode.HashTableCrefToSimVar crefToSimVarHT; - String errstr; - list subs; - Integer index; - case (cref, SimCode.SIMCODE(crefToSimVarHT = crefToSimVarHT) ) - algorithm - if BaseHashTable.hasKey(cref, crefToSimVarHT) then - sv := BaseHashTable.get(cref, crefToSimVarHT); - else - // lookup array variable and add offset for array element - if Flags.isSet(Flags.NF_SCALARIZE) then - sv := BaseHashTable.get(ComponentReference.crefStripLastSubs(cref), crefToSimVarHT); - subs := ComponentReference.crefLastSubs(cref); - sv.name := ComponentReference.crefSetLastSubs(sv.name, subs); - else - sv := BaseHashTable.get(ComponentReference.crefStripSubs(cref), crefToSimVarHT); - subs := ComponentReference.crefSubs(cref); - sv.name := ComponentReference.crefApplySubs(sv.name, subs); - end if; - sv.variable_index := match sv.variable_index - case SOME(index) - then SOME(index + getScalarElementIndex(subs, List.map(sv.numArrayElement, stringInt)) - 1); - end match; - end if; - sv := match sv.aliasvar - case SimCodeVar.NOALIAS() then sv; - case SimCodeVar.ALIAS(varName=cref) then cref2simvar(cref, simCode); /* Possibly not needed; can't really hurt that much though */ - case SimCodeVar.NEGATEDALIAS() then sv; - end match; - then sv; - - case (_, _) - equation - //print("cref2simvar: " + ComponentReference.printComponentRefStr(inCref) + " not found!\n"); - badcref = ComponentReference.makeCrefIdent("ERROR_cref2simvar_failed " + ComponentReference.printComponentRefStr(inCref), DAE.T_REAL_DEFAULT, {}); - then SimCodeVar.SIMVAR(badcref, BackendDAE.VARIABLE(), "", "", "", -2, NONE(), NONE(), NONE(), NONE(), false, DAE.T_REAL_DEFAULT, false, NONE(), SimCodeVar.NOALIAS(), DAE.emptyElementSource, SimCodeVar.INTERNAL(), NONE(), {}, false, true, false, NONE(), NONE()); - end matchcontinue; + try + SimCode.SIMCODE(crefToSimVarHT = crefToSimVarHT) := simCode; + outSimVar := simVarFromHT(inCref, crefToSimVarHT); + else + //print("cref2simvar: " + ComponentReference.printComponentRefStr(inCref) + " not found!\n"); + badcref := ComponentReference.makeCrefIdent("ERROR_cref2simvar_failed " + ComponentReference.printComponentRefStr(inCref), DAE.T_REAL_DEFAULT, {}); + outSimVar := SimCodeVar.SIMVAR(badcref, BackendDAE.VARIABLE(), "", "", "", -2, NONE(), NONE(), NONE(), NONE(), false, DAE.T_REAL_DEFAULT, false, NONE(), SimCodeVar.NOALIAS(), DAE.emptyElementSource, SimCodeVar.INTERNAL(), NONE(), {}, false, true, false, NONE(), NONE()); + end try; end cref2simvar; +public function simVarFromHT +"Used by templates to find SIMVAR for given cref (to gain representaion index info mainly)." + input DAE.ComponentRef inCref; + input HashTableCrefSimVar.HashTable crefToSimVarHT; + output SimCodeVar.SimVar outSimVar; +protected + DAE.ComponentRef cref, badcref; + SimCodeVar.SimVar sv; + list subs; + Integer index; +algorithm + try + if BaseHashTable.hasKey(inCref, crefToSimVarHT) then + sv := BaseHashTable.get(inCref, crefToSimVarHT); + else + // lookup array variable and add offset for array element + if Flags.isSet(Flags.NF_SCALARIZE) then + sv := BaseHashTable.get(ComponentReference.crefStripLastSubs(inCref), crefToSimVarHT); + subs := ComponentReference.crefLastSubs(inCref); + sv.name := ComponentReference.crefSetLastSubs(sv.name, subs); + else + sv := BaseHashTable.get(ComponentReference.crefStripSubs(inCref), crefToSimVarHT); + subs := ComponentReference.crefSubs(inCref); + sv.name := ComponentReference.crefApplySubs(sv.name, subs); + end if; + + sv.variable_index := match sv.variable_index + case SOME(index) + then SOME(index + getScalarElementIndex(subs, List.map(sv.numArrayElement, stringInt)) - 1); + end match; + end if; + sv := match sv.aliasvar + case SimCodeVar.NOALIAS() then sv; + case SimCodeVar.ALIAS(varName=cref) then simVarFromHT(cref, crefToSimVarHT); /* Possibly not needed; can't really hurt that much though */ + case SimCodeVar.NEGATEDALIAS() then sv; + end match; + else + //print("cref2simvar: " + ComponentReference.printComponentRefStr(inCref) + " not found!\n"); + badcref := ComponentReference.makeCrefIdent("ERROR_cref2simvar_failed " + ComponentReference.printComponentRefStr(inCref), DAE.T_REAL_DEFAULT, {}); + sv := SimCodeVar.SIMVAR(badcref, BackendDAE.VARIABLE(), "", "", "", -2, NONE(), NONE(), NONE(), NONE(), false, DAE.T_REAL_DEFAULT, false, NONE(), SimCodeVar.NOALIAS(), DAE.emptyElementSource, SimCodeVar.INTERNAL(), NONE(), {}, false, true, false, NONE(), NONE()); + end try; + outSimVar := sv; +end simVarFromHT; + +public function createJacContext + input Option jacHT; + output SimCodeFunction.Context outContext; +algorithm + outContext := SimCodeFunction.JACOBIAN_CONTEXT(jacHT); +end createJacContext; + public function codegenExpSanityCheck "Handle some things that Susan cannot handle: * Expand simulation context arrays that contain variables stored in different locations... * We could move collapsing arrays here since it should be safer to do so when we can lookup which index a variable corresponds to... diff --git a/Compiler/Stubs/SimCodeUtil.mo b/Compiler/Stubs/SimCodeUtil.mo index c45ca9094ac..6e4c0bc4a61 100644 --- a/Compiler/Stubs/SimCodeUtil.mo +++ b/Compiler/Stubs/SimCodeUtil.mo @@ -32,6 +32,14 @@ algorithm assert(false, getInstanceName()); end cref2simvar; +function simVarFromHT + input A inCref; + input B simCode; + output SimCodeVar.SimVar outSimVar; +algorithm + assert(false, getInstanceName()); +end simVarFromHT; + function codegenExpSanityCheck input output DAE.Exp e; input SimCodeFunction.Context context; diff --git a/Compiler/Template/CodegenC.tpl b/Compiler/Template/CodegenC.tpl index 1f4fdc8882c..011336f919e 100644 --- a/Compiler/Template/CodegenC.tpl +++ b/Compiler/Template/CodegenC.tpl @@ -160,16 +160,16 @@ end translateModel; extern const char* <%symbolName(modelNamePrefixStr,"zeroCrossingDescription")%>(int i, int **out_EquationIndexes); extern const char* <%symbolName(modelNamePrefixStr,"relationDescription")%>(int i); extern void <%symbolName(modelNamePrefixStr,"function_initSample")%>(DATA *data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianG")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianA")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianB")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianC")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianD")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"functionJacG_column")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"functionJacA_column")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"functionJacB_column")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"functionJacC_column")%>(void* data, threadData_t *threadData); - extern int <%symbolName(modelNamePrefixStr,"functionJacD_column")%>(void* data, threadData_t *threadData); + extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianG")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian); + extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianA")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian); + extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianB")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian); + extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianC")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian); + extern int <%symbolName(modelNamePrefixStr,"initialAnalyticJacobianD")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian); + extern int <%symbolName(modelNamePrefixStr,"functionJacG_column")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *thisJacobian, ANALYTIC_JACOBIAN *parentJacobian); + extern int <%symbolName(modelNamePrefixStr,"functionJacA_column")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *thisJacobian, ANALYTIC_JACOBIAN *parentJacobian); + extern int <%symbolName(modelNamePrefixStr,"functionJacB_column")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *thisJacobian, ANALYTIC_JACOBIAN *parentJacobian); + extern int <%symbolName(modelNamePrefixStr,"functionJacC_column")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *thisJacobian, ANALYTIC_JACOBIAN *parentJacobian); + extern int <%symbolName(modelNamePrefixStr,"functionJacD_column")%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *thisJacobian, ANALYTIC_JACOBIAN *parentJacobian); extern const char* <%symbolName(modelNamePrefixStr,"linear_model_frame")%>(void); extern const char* <%symbolName(modelNamePrefixStr,"linear_model_datarecovery_frame")%>(void); extern int <%symbolName(modelNamePrefixStr,"mayer")%>(DATA* data, modelica_real** res, short *); @@ -766,7 +766,7 @@ template simulationFile_jac_header(SimCode simCode) << /* Jacobians */ static const REAL_ATTRIBUTE dummyREAL_ATTRIBUTE = omc_dummyRealAttribute; - <%variableDefinitionsJacobians(jacobianMatrixes, modelNamePrefix(simCode))%> + <%symJacDefinition(jacobianMatrixes, modelNamePrefix(simCode))%> <%\n%> >> /* adrpo: leave a newline at the end of file to get rid of the warning */ @@ -1533,78 +1533,29 @@ template globalDataAliasVarArray(String _type, String _name, list items) end match end globalDataAliasVarArray; -template variableDefinitionsJacobians(list JacobianMatrixes, String modelNamePrefix) "template variableDefinitionsJacobians +template symJacDefinition(list JacobianMatrixes, String modelNamePrefix) "template variableDefinitionsJacobians Generates defines for jacobian vars." ::= - let analyticVars = (JacobianMatrixes |> JAC_MATRIX(columns=jacColumn, seedVars=seedVars, matrixName=name, jacobianIndex=indexJacobian) => - let varsDef = variableDefinitionsJacobians2(indexJacobian, jacColumn, seedVars, name) + let symbolicJacsDefine = (JacobianMatrixes |> JAC_MATRIX(columns=jacColumn, seedVars=seedVars, matrixName=name, jacobianIndex=indexJacobian) => << #if defined(__cplusplus) extern "C" { #endif #define <%symbolName(modelNamePrefix,"INDEX_JAC_")%><%name%> <%indexJacobian%> - int <%symbolName(modelNamePrefix,"functionJac")%><%name%>_column(void* data, threadData_t *threadData); - int <%symbolName(modelNamePrefix,"initialAnalyticJacobian")%><%name%>(void* data, threadData_t *threadData); + int <%symbolName(modelNamePrefix,"functionJac")%><%name%>_column(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *thisJacobian, ANALYTIC_JACOBIAN *parentJacobian); + int <%symbolName(modelNamePrefix,"initialAnalyticJacobian")%><%name%>(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian); #if defined(__cplusplus) } #endif - <%varsDef%> >> ;separator="\n";empty) << /* Jacobian Variables */ - <%analyticVars%> + <%symbolicJacsDefine%> >> -end variableDefinitionsJacobians; - -template variableDefinitionsJacobians2(Integer indexJacobian, list jacobianColumn, list seedVars, String name) "template variableDefinitionsJacobians2 - Generates Matrixes for Linear Model." -::= - let seedVarsResult = (seedVars |> var hasindex index0 => - jacobianVarDefine(var, indexJacobian, index0, name) - ;separator="\n") - let columnVarsResult = (jacobianColumn |> JAC_COLUMN(columnVars=vars) => - (vars |> var hasindex index0 => jacobianVarDefine(var, indexJacobian, index0, name);separator="\n") - ;separator="\n\n") - /* generate at least one print command to have the same index and avoid the strange side effect */ - << - /* <%name%> */ - <%seedVarsResult%> - <%columnVarsResult%> - >> -end variableDefinitionsJacobians2; - -template jacobianVarDefine(SimVar simVar, Integer indexJac, Integer index0, String matrixName) "template jacobianVarDefine - " -::= - match simVar - case SIMVAR(arrayCref=arrayCref,aliasvar=NOALIAS(),name=name,index=index) then - let crefName = crefDefine(name) - let typeName = match varKind - case BackendDAE.JAC_VAR() then 'resultVars[<%index%>]' - case BackendDAE.JAC_DIFF_VAR() then 'tmpVars[<%index%>]' - case BackendDAE.SEED_VAR() then 'seedVars[<%index0%>]' - else 'UNKNOWN KIND' - let arrayCrefString = if Util.isSome(arrayCref) then cref(Util.getOption(arrayCref)) else '' - let arrayDefine = if Util.isSome(arrayCref) then - << - /* array <%arrayCrefString%> */ - #define _<%arrayCrefString%>(i) data->simulationInfo->analyticJacobians[<%indexJac%>].<%typeName%> - #define <%arrayCrefString%> _<%arrayCrefString%>(0) - >> else '' - << - <%arrayDefine%> - /* <%crefName%> */ - #define _<%crefName%>(i) data->simulationInfo->analyticJacobians[<%indexJac%>].<%typeName%> - #define <%crefName%> _<%crefName%>(0) - <%if stringEq('<%crefName%>', '$P<%BackendDAE.optimizationMayerTermName%>$P$pDERC$PdummyVarC') then "\n"+'#define <%crefName%>$indexdiffed <%index%>' - %><%if stringEq('<%crefName%>', '$P<%BackendDAE.optimizationLagrangeTermName%>$P$pDERB$PdummyVarB') then "\n"+'#define <%crefName%>$indexdiffed <%index%>' - %><%if stringEq('<%crefName%>', '$P<%BackendDAE.optimizationLagrangeTermName%>$P$pDERC$PdummyVarC') then "\n"+'#define <%crefName%>$indexdiffed <%index%>'%> - >> - end match -end jacobianVarDefine; +end symJacDefinition; template aliasVarNameType(AliasVariable var) "Generates type of alias." @@ -1939,6 +1890,7 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].setA = setLinearMatrixA<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].setb = setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; >> case SOME(__) then let size = listLength(ls.vars) @@ -1960,6 +1912,7 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; >> else error(sourceInfo(), ' No jacobian create for linear system <%ls.index%>.') @@ -1986,6 +1939,7 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].setA = setLinearMatrixA<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].setb = setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; + linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; assertStreamPrint(NULL, nLinearSystems > <%at.indexLinearSystem%>, "Internal Error: nLinearSystems mismatch!"); linearSystemData[<%at.indexLinearSystem%>].equationIndex = <%at.index%>; @@ -1996,6 +1950,7 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%at.indexLinearSystem%>].setA = setLinearMatrixA<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].setb = setLinearVectorb<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%at.index%>; + linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; linearSystemData[<%at.indexLinearSystem%>].checkConstraints = checkConstraints<%at.index%>; >> case SOME(__) then @@ -2021,6 +1976,7 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%ls.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac%>; linearSystemData[<%ls.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac%>; linearSystemData[<%ls.indexLinearSystem%>].jacobianIndex = <%jacIndex%>; + linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; linearSystemData[<%ls.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%ls.index%>; linearSystemData[<%ls.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%ls.index%>; @@ -2035,6 +1991,7 @@ template functionInitialLinearSystemsTemp(list linearSystems, Strin linearSystemData[<%at.indexLinearSystem%>].analyticalJacobianColumn = <%generatedJac2%>; linearSystemData[<%at.indexLinearSystem%>].initialAnalyticalJacobian = <%initialJac2%>; linearSystemData[<%at.indexLinearSystem%>].jacobianIndex = <%jacIndex2%>; + linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = NULL; linearSystemData[<%at.indexLinearSystem%>].setA = NULL;//setLinearMatrixA<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].setb = NULL; //setLinearVectorb<%at.index%>; linearSystemData[<%at.indexLinearSystem%>].initializeStaticLSData = initializeStaticLSData<%at.index%>; @@ -2070,12 +2027,19 @@ template functionSetupLinearSystemsTemp(list linearSystems, String let &auxFunction = buffer "" let &tmp = buffer "" let xlocs = (ls.vars |> var hasindex i0 => '<%cref(varName(var))%> = xloc[<%i0%>];' ;separator="\n") - let prebody = (ls.residual |> eq2 => + let prebody = (match ls.partOfJac + case true then + (ls.residual |> eq2 => + functionExtraResidualsPreBodyJacobian(eq2, &tmp, modelNamePrefix) + ;separator="\n") + case false then + (ls.residual |> eq2 => functionExtraResidualsPreBody(eq2, &tmp, modelNamePrefix) - ;separator="\n") + ;separator="\n") + end match) let body = (ls.residual |> eq2 as SES_RESIDUAL(__) hasindex i0 => - let &preExp = buffer "" /*BUFD*/ - let expPart = daeExp(eq2.exp, contextSimulationDiscrete, + let &preExp = buffer "" /*BUFD*/ + let expPart = daeExp(eq2.exp, contextSimulationDiscrete, &preExp /*BUFC*/, &varDeclsRes, &auxFunction) << <% if profileAll() then 'SIM_PROF_TICK_EQ(<%eq2.index%>);' %> @@ -2099,6 +2063,8 @@ template functionSetupLinearSystemsTemp(list linearSystems, String DATA *data = (DATA*) ((void**)dataIn[0]); threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]); const int equationIndexes[2] = {1,<%ls.index%>}; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%> + ANALYTIC_JACOBIAN* jacobian = NULL; <%varDeclsRes%> <% if profileAll() then 'SIM_PROF_TICK_EQ(<%ls.index%>);' %> <% if profileSome() then 'SIM_PROF_ADD_NCALL_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex,1);' %> @@ -2147,6 +2113,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> <%varDecls%> <%MatrixA%> } @@ -2156,6 +2123,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> <%varDecls2%> <%vectorb%> } @@ -2233,6 +2201,8 @@ template functionSetupLinearSystemsTemp(list linearSystems, String DATA *data = (DATA*) ((void**)dataIn[0]); threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]); const int equationIndexes[2] = {1,<%ls.index%>}; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%> + ANALYTIC_JACOBIAN* jacobian = NULL; <%varDeclsRes%> <% if profileAll() then 'SIM_PROF_TICK_EQ(<%ls.index%>);' %> <% if profileSome() then 'SIM_PROF_ADD_NCALL_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex,1);' %> @@ -2260,6 +2230,8 @@ template functionSetupLinearSystemsTemp(list linearSystems, String DATA *data = (DATA*) ((void**)dataIn[0]); threadData_t *threadData = (threadData_t*) ((void**)dataIn[1]); const int equationIndexes[2] = {1,<%at.index%>}; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian;'%> + ANALYTIC_JACOBIAN* jacobian = NULL; <%varDeclsRes2%> <% if profileAll() then 'SIM_PROF_TICK_EQ(<%at.index%>);' %> <% if profileSome() then 'SIM_PROF_ADD_NCALL_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%at.index%>).profileBlockIndex,1);' %> @@ -2332,6 +2304,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> <%varDecls%> <%MatrixA%> } @@ -2341,6 +2314,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%ls.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> <%varDecls2%> <%vectorb%> } @@ -2360,6 +2334,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%at.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> <%varDecls3%> <%MatrixA2%> } @@ -2369,6 +2344,7 @@ template functionSetupLinearSystemsTemp(list linearSystems, String const int equationIndexes[2] = {1,<%at.index%>}; DATA* data = (DATA*) inData; LINEAR_SYSTEM_DATA* linearSystemData = (LINEAR_SYSTEM_DATA*) systemData; + <% if ls.partOfJac then 'ANALYTIC_JACOBIAN* parentJacobian = linearSystemData->parentJacobian;'%> <%varDecls4%> <%vectorb2%> } @@ -2539,6 +2515,23 @@ template functionExtraResidualsPreBody(SimEqSystem eq, Text &eqs, String modelNa end match end functionExtraResidualsPreBody; +template functionExtraResidualsPreBodyJacobian(SimEqSystem eq, Text &eqs, String modelNamePrefixStr) + "Generates an equation." +::= + let () = tmpTickReset(0) + match eq + case e as SES_RESIDUAL(__) + then "" + else + let &eqs += equation_impl(-1, eq, contextJacobian, modelNamePrefixStr) + << + /* local constraints */ + <%createLocalConstraints(eq)%> + <%equation_callJacobian(eq, modelNamePrefixStr)%> + >> + end match +end functionExtraResidualsPreBodyJacobian; + template createLocalConstraints(SimEqSystem eq) "Generates an equation." ::= @@ -4642,8 +4635,8 @@ template functionAnalyticJacobians(list JacobianMatrixes,String ::= let initialjacMats = (JacobianMatrixes |> JAC_MATRIX(columns=mat, seedVars=vars, matrixName=name, sparsity=sparsepattern, coloredCols=colorList, maxColorCols=maxColor, jacobianIndex=indexJacobian) => initialAnalyticJacobians(mat, vars, name, sparsepattern, colorList, maxColor, modelNamePrefix); separator="\n") - let jacMats = (JacobianMatrixes |> JAC_MATRIX(columns=mat, seedVars=vars, matrixName=name, partitionIndex=partIdx) => - generateMatrix(mat, vars, name, partIdx, modelNamePrefix) ;separator="\n") + let jacMats = (JacobianMatrixes |> JAC_MATRIX(columns=mat, seedVars=vars, matrixName=name, partitionIndex=partIdx, crefsHT=crefsHT) => + generateMatrix(mat, vars, name, partIdx, crefsHT, modelNamePrefix) ;separator="\n") << <%initialjacMats%> @@ -4660,7 +4653,7 @@ template initialAnalyticJacobians(list jacobianColumn, list<%matrixname%>(void* inData, threadData_t *threadData) + int <%symbolName(modelNamePrefix,"initialAnalyticJacobian")%><%matrixname%>(void* inData, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian) { TRACE_PUSH TRACE_POP @@ -4673,42 +4666,40 @@ match sparsepattern let sizeleadindex = listLength(sparsepattern) let colPtr = genSPCRSPtr(listLength(sparsepattern), sparsepattern, "colPtrIndex") let rowIndex = genSPCRSRows(lengthListElements(unzipSecond(sparsepattern)), sparsepattern, "rowIndex") - let colorString = genSPColors(colorList, "data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols") + let colorString = genSPColors(colorList, "jacobian->sparsePattern.colorCols") let indexColumn = (jacobianColumn |> JAC_COLUMN(numberOfResultVars=n) => '<%n%>';separator="\n") let tmpvarsSize = (jacobianColumn |> JAC_COLUMN(columnVars=vars) => listLength(vars);separator="\n") let index_ = listLength(seedVars) << OMC_DISABLE_OPT - int <%symbolName(modelNamePrefix,"initialAnalyticJacobian")%><%matrixname%>(void* inData, threadData_t *threadData) + int <%symbolName(modelNamePrefix,"initialAnalyticJacobian")%><%matrixname%>(void* inData, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian) { TRACE_PUSH DATA* data = ((DATA*)inData); - int index = <%symbolName(modelNamePrefix,"INDEX_JAC_")%><%matrixname%>; <%colPtr%> <%rowIndex%> int i = 0; - data->simulationInfo->analyticJacobians[index].sizeCols = <%index_%>; - data->simulationInfo->analyticJacobians[index].sizeRows = <%indexColumn%>; - data->simulationInfo->analyticJacobians[index].sizeTmpVars = <%tmpvarsSize%>; - data->simulationInfo->analyticJacobians[index].seedVars = (modelica_real*) calloc(<%index_%>,sizeof(modelica_real)); - data->simulationInfo->analyticJacobians[index].resultVars = (modelica_real*) calloc(<%indexColumn%>,sizeof(modelica_real)); - data->simulationInfo->analyticJacobians[index].tmpVars = (modelica_real*) calloc(<%tmpvarsSize%>,sizeof(modelica_real)); - data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex = (unsigned int*) malloc((<%sizeleadindex%>+1)*sizeof(int)); - data->simulationInfo->analyticJacobians[index].sparsePattern.index = (unsigned int*) malloc(<%sp_size_index%>*sizeof(int)); - data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros = <%sp_size_index%>; - data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols = (unsigned int*) malloc(<%index_%>*sizeof(int)); - data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors = <%maxColor%>; - data->simulationInfo->analyticJacobians[index].jacobian = NULL; + jacobian->sizeCols = <%index_%>; + jacobian->sizeRows = <%indexColumn%>; + jacobian->sizeTmpVars = <%tmpvarsSize%>; + jacobian->seedVars = (modelica_real*) calloc(<%index_%>,sizeof(modelica_real)); + jacobian->resultVars = (modelica_real*) calloc(<%indexColumn%>,sizeof(modelica_real)); + jacobian->tmpVars = (modelica_real*) calloc(<%tmpvarsSize%>,sizeof(modelica_real)); + jacobian->sparsePattern.leadindex = (unsigned int*) malloc((<%sizeleadindex%>+1)*sizeof(int)); + jacobian->sparsePattern.index = (unsigned int*) malloc(<%sp_size_index%>*sizeof(int)); + jacobian->sparsePattern.numberOfNoneZeros = <%sp_size_index%>; + jacobian->sparsePattern.colorCols = (unsigned int*) malloc(<%index_%>*sizeof(int)); + jacobian->sparsePattern.maxColors = <%maxColor%>; /* write lead index of compressed sparse column */ - memcpy(data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex, colPtrIndex, (<%sizeleadindex%>+1)*sizeof(int)); + memcpy(jacobian->sparsePattern.leadindex, colPtrIndex, (<%sizeleadindex%>+1)*sizeof(int)); for(i=2;i<<%sizeleadindex%>+1;++i) - data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i] += data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i-1]; + jacobian->sparsePattern.leadindex[i] += jacobian->sparsePattern.leadindex[i-1]; /* call sparse index */ - memcpy(data->simulationInfo->analyticJacobians[index].sparsePattern.index, rowIndex, <%sp_size_index%>*sizeof(int)); + memcpy(jacobian->sparsePattern.index, rowIndex, <%sp_size_index%>*sizeof(int)); /* write color array */ <%colorString%> @@ -4719,7 +4710,7 @@ match sparsepattern end match end initialAnalyticJacobians; -template generateMatrix(list jacobianColumn, list seedVars, String matrixname, Integer partIdx, String modelNamePrefix) +template generateMatrix(list jacobianColumn, list seedVars, String matrixname, Integer partIdx, Option jacHT, String modelNamePrefix) "This template generates source code for a single jacobian in dense format and sparse format. This is a helper of template functionAnalyticJacobians" ::= @@ -4727,7 +4718,7 @@ template generateMatrix(list jacobianColumn, list seedVa match nRows case "0" then << - int <%symbolName(modelNamePrefix,"functionJac")%><%matrixname%>_column(void* data, threadData_t *threadData) + int <%symbolName(modelNamePrefix,"functionJac")%><%matrixname%>_column(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian, ANALYTIC_JACOBIAN *parentJacobian) { TRACE_PUSH TRACE_POP @@ -4738,7 +4729,7 @@ template generateMatrix(list jacobianColumn, list seedVa match seedVars case {} then << - int <%symbolName(modelNamePrefix,"functionJac")%><%matrixname%>_column(void* data, threadData_t *threadData) + int <%symbolName(modelNamePrefix,"functionJac")%><%matrixname%>_column(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian, ANALYTIC_JACOBIAN *parentJacobian) { TRACE_PUSH TRACE_POP @@ -4747,7 +4738,7 @@ template generateMatrix(list jacobianColumn, list seedVa >> case _ then let jacMats = (jacobianColumn |> JAC_COLUMN(columnEqns=eqs) => - functionJac(eqs, partIdx, matrixname, modelNamePrefix) + functionJac(eqs, partIdx, matrixname, jacHT, modelNamePrefix) ;separator="\n") let indexColumn = (jacobianColumn |> JAC_COLUMN(numberOfResultVars=nRows) => nRows @@ -4759,21 +4750,22 @@ template generateMatrix(list jacobianColumn, list seedVa end match end generateMatrix; -template functionJac(list jacEquations, Integer partIdx, String matrixName, String modelNamePrefix) "template functionJac +template functionJac(list jacEquations, Integer partIdx, String matrixName, Option jacHT, String modelNamePrefix) "template functionJac This template generates functions for each column of a single jacobian. This is a helper of generateMatrix." ::= + << <%(jacEquations |> eq => - equation_impl(partIdx, eq, contextSimulationNonDiscrete, modelNamePrefix); separator="\n")%> - int <%symbolName(modelNamePrefix,"functionJac")%><%matrixName%>_column(void* inData, threadData_t *threadData) + equation_impl(partIdx, eq, createJacContext(jacHT), modelNamePrefix); separator="\n")%> + int <%symbolName(modelNamePrefix,"functionJac")%><%matrixName%>_column(void* inData, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian, ANALYTIC_JACOBIAN *parentJacobian) { TRACE_PUSH DATA* data = ((DATA*)inData); int index = <%symbolName(modelNamePrefix,"INDEX_JAC_")%><%matrixName%>; <%(jacEquations |> eq => - equation_call(eq, modelNamePrefix); separator="\n")%> + equation_callJacobian(eq, modelNamePrefix); separator="\n")%> TRACE_POP return 0; @@ -4925,120 +4917,152 @@ template equation_impl2(Integer clockIndex, SimEqSystem eq, Context context, Str ::= let OMC_NO_OPT = if noOpt then 'OMC_DISABLE_OPT<%\n%>' else (match eq case SES_LINEAR(__) then 'OMC_DISABLE_OPT<%\n%>') match eq - case SES_ALIAS(__) then 'extern void <%symbolName(modelNamePrefix,"eqFunction")%>_<%aliasOf%>(DATA *data, threadData_t *threadData);<%\n%>' - case SES_ALGORITHM(statements={}) - then "" - else - ( - let ix = equationIndex(eq) /*System.tmpTickIndex(10)*/ - let ix2 = match eq - case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) - case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__), alternativeTearing = SOME(at as NONLINEARSYSTEM(__))) then - equationIndexAlternativeTearing(eq) - else - "" - let &tmp = buffer "" - let &varD = buffer "" - let &tempeqns = buffer "" - let &tempeqns2 = buffer "" - let() = System.tmpTickResetIndex(0,1) /* Boxed array indices */ - let disc = match context - case SIMULATION_CONTEXT(genDiscrete=true) then 1 - else 0 - let x = match eq - case e as SES_SIMPLE_ASSIGN(__) - case e as SES_SIMPLE_ASSIGN_CONSTRAINTS(__) - then equationSimpleAssign(e, context, &varD, &tempeqns) - case e as SES_ARRAY_CALL_ASSIGN(__) - then equationArrayCallAssign(e, context, &varD, &tempeqns) - case e as SES_IFEQUATION(__) - then equationIfEquationAssign(e, context, &varD, &tempeqns, modelNamePrefix) - case e as SES_INVERSE_ALGORITHM(__) - case e as SES_ALGORITHM(__) - then equationAlgorithm(e, context, &varD, &tempeqns) - case e as SES_LINEAR(__) - then equationLinear(e, context, &varD) - case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__)) then - let &tempeqns += (nls.eqs |> eq => 'void <%symbolName(modelNamePrefix,"eqFunction")%>_<%equationIndex(eq)%>(DATA*,threadData_t*);' ; separator = "\n") - equationNonlinear(e, context, modelNamePrefix) - case e as SES_WHEN(__) - then equationWhen(e, context, &varD, &tempeqns) - case e as SES_RESIDUAL(__) - then "NOT IMPLEMENTED EQUATION SES_RESIDUAL" - case e as SES_MIXED(__) - then - let &eqs = buffer "" - let res = equationMixed(e, context, &eqs, modelNamePrefix) - eqs + res - case e as SES_FOR_LOOP(__) - then equationForLoop(e, context, &varD, &tempeqns) - else - error(sourceInfo(), "NOT IMPLEMENTED EQUATION equation_") - let x2 = match eq - case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) then - equationLinearAlternativeTearing(e, context, &varD) - case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__), alternativeTearing = SOME(at as NONLINEARSYSTEM(__))) then - let &tempeqns2 += (at.eqs |> eq => 'void <%symbolName(modelNamePrefix,"eqFunction")%>_<%equationIndex(eq)%>(DATA*,threadData_t*);' ; separator = "\n") - equationNonlinearAlternativeTearing(e, context, modelNamePrefix) - else - "" - let &varD += addRootsTempArray() - let clockIndex_ = if intLt(clockIndex, 0) then '' else 'const int clockIndex = <%clockIndex%>;' + case SES_ALIAS(__) then 'extern void <%symbolName(modelNamePrefix,"eqFunction")%>_<%aliasOf%>(DATA *data, threadData_t *threadData);<%\n%>' + case e as SES_ALGORITHM(statements={}) then "" + else + ( + let ix = equationIndex(eq) /*System.tmpTickIndex(10)*/ + let ix2 = match eq + case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) + case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__), alternativeTearing = SOME(at as NONLINEARSYSTEM(__))) + then equationIndexAlternativeTearing(eq) + else "" + let &tmp = buffer "" + let &varD = buffer "" + let &tempeqns = buffer "" + let &tempeqns2 = buffer "" + let() = System.tmpTickResetIndex(0,1) /* Boxed array indices */ + let disc = match context case SIMULATION_CONTEXT(genDiscrete=true) then 1 else 0 + let x = match eq + + case e as SES_SIMPLE_ASSIGN(__) + case e as SES_SIMPLE_ASSIGN_CONSTRAINTS(__) + then equationSimpleAssign(e, context, &varD, &tempeqns) + + case e as SES_ARRAY_CALL_ASSIGN(__) + then equationArrayCallAssign(e, context, &varD, &tempeqns) + + case e as SES_IFEQUATION(__) + then equationIfEquationAssign(e, context, &varD, &tempeqns, modelNamePrefix) + + case e as SES_ALGORITHM(__) + then equationAlgorithm(e, context, &varD, &tempeqns) + + case e as SES_INVERSE_ALGORITHM(__) + then equationAlgorithm(e, context, &varD, &tempeqns) + + case e as SES_LINEAR(__) + then equationLinear(e, context, &varD) + + case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__)) + then + let &tempeqns += (nls.eqs |> eq => 'void <%symbolName(modelNamePrefix,"eqFunction")%>_<%equationIndex(eq)%>(DATA*, threadData_t*);' ; separator = "\n") + equationNonlinear(e, context, modelNamePrefix) + + case e as SES_WHEN(__) + then equationWhen(e, context, &varD, &tempeqns) + + case e as SES_RESIDUAL(__) + then "NOT IMPLEMENTED EQUATION SES_RESIDUAL" + + case e as SES_MIXED(__) + then + let &eqs = buffer "" + let res = equationMixed(e, context, &eqs, modelNamePrefix) + eqs + res + + case e as SES_FOR_LOOP(__) + then equationForLoop(e, context, &varD, &tempeqns) + + else "NOT IMPLEMENTED EQUATION equation_" + + let x2 = match eq + case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) + then equationLinearAlternativeTearing(e, context, &varD) + + case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__), alternativeTearing = SOME(at as NONLINEARSYSTEM(__))) + then + let &tempeqns2 += (at.eqs |> eq => 'void <%symbolName(modelNamePrefix,"eqFunction")%>_<%equationIndex(eq)%>(DATA*,threadData_t*);' ; separator = "\n") + equationNonlinearAlternativeTearing(e, context, modelNamePrefix) + else "" + + let &varD += addRootsTempArray() + let clockIndex_ = if intLt(clockIndex, 0) then '' else 'const int clockIndex = <%clockIndex%>;' + + match eq + // dynamic tearing + case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) + case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__), alternativeTearing = SOME(at as NONLINEARSYSTEM(__))) then + << - match eq - // dynamic tearing - case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) - case e as SES_NONLINEAR(nlSystem=nls as NONLINEARSYSTEM(__), alternativeTearing = SOME(at as NONLINEARSYSTEM(__))) then - << + <%tempeqns%> + /* + <%dumpEqs(fill(eq,1))%> + */ + <%OMC_NO_OPT%><% if static then "static "%>int <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix%>(DATA *data, threadData_t *threadData) + { + TRACE_PUSH + <%clockIndex_%> + const int equationIndexes[2] = {1,<%ix%>}; + <%&varD%> + <%x%> + TRACE_POP + } - <%tempeqns%> - /* - <%dumpEqs(fill(eq,1))%> - */ - <%OMC_NO_OPT%><% if static then "static "%>int <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix%>(DATA *data, threadData_t *threadData) - { - TRACE_PUSH - <%clockIndex_%> - const int equationIndexes[2] = {1,<%ix%>}; - <%&varD%> - <%x%> - TRACE_POP - } + <%tempeqns2%> + /* + <%dumpEqsAlternativeTearing(fill(eq,1))%> + */ + <%OMC_NO_OPT%><% if static then "static "%>void <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix2%>(DATA *data, threadData_t *threadData) + { + TRACE_PUSH + <%clockIndex_%> + const int equationIndexes[2] = {1,<%ix2%>}; + <%&varD%> + <%x2%> + TRACE_POP + } + >> - <%tempeqns2%> - /* - <%dumpEqsAlternativeTearing(fill(eq,1))%> - */ - <%OMC_NO_OPT%><% if static then "static "%>void <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix2%>(DATA *data, threadData_t *threadData) - { - TRACE_PUSH - <%clockIndex_%> - const int equationIndexes[2] = {1,<%ix2%>}; - <%&varD%> - <%x2%> - TRACE_POP - } - >> + // no dynamic tearing + else + match context + case JACOBIAN_CONTEXT() + then + << - // no dynamic tearing - else - << + <%tempeqns%> + /* + <%dumpEqs(fill(eq,1))%> + */ + <%OMC_NO_OPT%><% if static then "static "%>void <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix%>(DATA *data, threadData_t *threadData, ANALYTIC_JACOBIAN *jacobian, ANALYTIC_JACOBIAN *parentJacobian) + { + TRACE_PUSH + <%clockIndex_%> + const int equationIndexes[2] = {1,<%ix%>}; + <%&varD%> + <%x%> + TRACE_POP + } + >> + else + << - <%tempeqns%> - /* - <%dumpEqs(fill(eq,1))%> - */ - <%OMC_NO_OPT%><% if static then "static "%>void <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix%>(DATA *data, threadData_t *threadData) - { - TRACE_PUSH - <%clockIndex_%> - const int equationIndexes[2] = {1,<%ix%>}; - <%&varD%> - <%x%> - TRACE_POP - } - >> + <%tempeqns%> + /* + <%dumpEqs(fill(eq,1))%> + */ + <%OMC_NO_OPT%><% if static then "static "%>void <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix%>(DATA *data, threadData_t *threadData) + { + TRACE_PUSH + <%clockIndex_%> + const int equationIndexes[2] = {1,<%ix%>}; + <%&varD%> + <%x%> + TRACE_POP + } + >> ) end equation_impl2; @@ -5067,6 +5091,30 @@ template equation_call(SimEqSystem eq, String modelNamePrefix) ) end equation_call; +template equation_callJacobian(SimEqSystem eq, String modelNamePrefix) + "Generates the equation calls for jacobians fucntion." +::= + match eq + case e as SES_ALGORITHM(statements={}) + then "" + else + ( + let ix = match eq + case e as SES_LINEAR(alternativeTearing = SOME(LINEARSYSTEM)) + case e as SES_NONLINEAR(alternativeTearing = SOME(NONLINEARSYSTEM)) then + equationIndexAlternativeTearing(eq) + else + equationIndex(eq) + end match + << + <% if profileAll() then 'SIM_PROF_TICK_EQ(<%ix%>);' %> + <%symbolName(modelNamePrefix,"eqFunction")%>_<%ix%>(data, threadData, jacobian, parentJacobian); + <% if profileAll() then 'SIM_PROF_ACC_EQ(<%ix%>);' %> + >> + ) +end equation_callJacobian; + + template equationForward_(SimEqSystem eq, Context context, String modelNamePrefixStr) "Generates an equation. This template should not be used for a SES_RESIDUAL. @@ -5135,7 +5183,7 @@ case SES_SIMPLE_ASSIGN_CONSTRAINTS(__) then << <%modelicaLine(eqInfo(eq))%> <%preExp%> - <%cref(cref)%> = <%expPart%>; + <%contextCref(cref, context, auxFunction)%> = <%expPart%>; <%postExp%> <%endModelicaLine()%> >> @@ -5244,17 +5292,19 @@ match eq case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = at) then let returnval = match at case at as SOME(__) then 'return 1;' case at as NONE() then '' let returnval2 = match at case at as SOME(__) then 'return 0;' case at as NONE() then '' + let auxFunctions = "" << /* Linear equation system */ int retValue; + double aux_x[<%listLength(ls.vars)%>] = { <%ls.vars |> SIMVAR(__) hasindex i0 => '<%contextCrefOld(name, context, auxFunctions, 1)%>' ;separator=","%> }; if(ACTIVE_STREAM(LOG_DT)) { infoStreamPrint(LOG_DT, 1, "Solving linear system <%ls.index%> (STRICT TEARING SET if tearing enabled) at time = %18.10e", data->localData[0]->timeValue); messageClose(LOG_DT); } <% if profileSome() then 'SIM_PROF_TICK_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex);' %> - <%ls.vars |> SIMVAR(__) hasindex i0 => 'data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].x[<%i0%>] = <%crefOld(name,1)%>;' ;separator="\n"%> - retValue = solve_linear_system(data, threadData, <%ls.indexLinearSystem%>); + <% if ls.partOfJac then 'data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].parentJacobian = jacobian;'%> + retValue = solve_linear_system(data, threadData, <%ls.indexLinearSystem%>, &aux_x[0]); /* check if solution process was successful */ if (retValue > 0){ @@ -5263,7 +5313,7 @@ case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = at) th <%returnval2%> } /* write solution */ - <%ls.vars |> SIMVAR(__) hasindex i0 => '<%cref(name)%> = data->simulationInfo->linearSystemData[<%ls.indexLinearSystem%>].x[<%i0%>];' ;separator="\n"%> + <%ls.vars |> SIMVAR(__) hasindex i0 => '<%contextCref(name, context, auxFunctions)%> = aux_x[<%i0%>];' ;separator="\n"%> <% if profileSome() then 'SIM_PROF_ACC_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%ls.index%>).profileBlockIndex);' %> <%returnval%> >> @@ -5275,9 +5325,11 @@ template equationLinearAlternativeTearing(SimEqSystem eq, Context context, Text ::= match eq case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(at as LINEARSYSTEM(__))) then + let auxFunctions = "" << /* Linear equation system */ int retValue; + if(ACTIVE_STREAM(LOG_DT)) { infoStreamPrint(LOG_DT, 1, "Solving linear system <%at.index%> (CASUAL TEARING SET, strict: <%ls.index%>) at time = %18.10e", data->localData[0]->timeValue); @@ -5286,12 +5338,12 @@ case e as SES_LINEAR(lSystem=ls as LINEARSYSTEM(__), alternativeTearing = SOME(a <% if profileSome() then 'SIM_PROF_TICK_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%at.index%>).profileBlockIndex);' %> if (data->simulationInfo->linearSystemData[<%at.indexLinearSystem%>].checkConstraints(data, threadData) == 1) { - <%at.vars |> SIMVAR(__) hasindex i0 => 'data->simulationInfo->linearSystemData[<%at.indexLinearSystem%>].x[<%i0%>] = <%crefOld(name,1)%>;' ;separator="\n"%> - retValue = solve_linear_system(data, threadData, <%at.indexLinearSystem%>); + double aux_x[<%listLength(at.vars)%>] = { <%at.vars |> SIMVAR(__) hasindex i0 => '<%contextCref(name, context, auxFunctions)%>' ;separator=","%> }; + retValue = solve_linear_system(data, threadData, <%at.indexLinearSystem%>, &aux_x[0]); /* The casual tearing set found a solution */ if (retValue == 0){ /* write solution */ - <%at.vars |> SIMVAR(__) hasindex i0 => '<%cref(name)%> = data->simulationInfo->linearSystemData[<%at.indexLinearSystem%>].x[<%i0%>];' ;separator="\n"%> + <%at.vars |> SIMVAR(__) hasindex i0 => '<%contextCref(name, context, auxFunctions)%> = aux_x[<%i0%>];' ;separator="\n"%> <% if profileSome() then 'SIM_PROF_ACC_EQ(modelInfoGetEquation(&data->modelData->modelDataXml,<%at.index%>).profileBlockIndex);' %> } } @@ -5967,12 +6019,14 @@ template optimizationComponents1(ClassAttributes classAttribute, SimCode simCode let objectiveFunction = match objetiveE case SOME(exp) then + let setMayerIndex = if stringEq(crefToIndex(createDifferentiatedCrefName(makeUntypedCrefIdent(BackendDAE.optimizationMayerTermName), makeUntypedCrefIdent("dummyVarC"), "C")), "-2") then '' else << - *res = &<%cref(makeUntypedCrefIdent(BackendDAE.optimizationMayerTermName))%>; - #ifdef $P<%BackendDAE.optimizationMayerTermName%>$P$pDERC$PdummyVarC$indexdiffed - *index_Dres = $P<%BackendDAE.optimizationMayerTermName%>$P$pDERC$PdummyVarC$indexdiffed; + *index_Dres = <%crefToIndex(createDifferentiatedCrefName(makeUntypedCrefIdent(BackendDAE.optimizationMayerTermName), makeUntypedCrefIdent("dummyVarC"), "C"))%>; return 0; - #endif + >> + << + *res = &<%cref(makeUntypedCrefIdent(BackendDAE.optimizationMayerTermName))%>; + <%setMayerIndex%> >> let startTimeOpt = match startTimeE @@ -5984,13 +6038,15 @@ template optimizationComponents1(ClassAttributes classAttribute, SimCode simCode let objectiveIntegrand = match objectiveIntegrandE case SOME(exp) then + let setLagrangeIndex = if stringEq(crefToIndex(createDifferentiatedCrefName(makeUntypedCrefIdent(BackendDAE.optimizationLagrangeTermName), makeUntypedCrefIdent("dummyVarB"), "B")), "-2") then '' else << - *res = &<%cref(makeUntypedCrefIdent(BackendDAE.optimizationLagrangeTermName))%>; - #ifdef $P<%BackendDAE.optimizationLagrangeTermName%>$P$pDERB$PdummyVarB$indexdiffed - *index_DresB = $P<%BackendDAE.optimizationLagrangeTermName%>$P$pDERB$PdummyVarB$indexdiffed; - *index_DresC = $P<%BackendDAE.optimizationLagrangeTermName%>$P$pDERC$PdummyVarC$indexdiffed; + *index_DresB = <%crefToIndex(createDifferentiatedCrefName(makeUntypedCrefIdent(BackendDAE.optimizationLagrangeTermName), makeUntypedCrefIdent("dummyVarB"), "B"))%>; + *index_DresC = <%crefToIndex(createDifferentiatedCrefName(makeUntypedCrefIdent(BackendDAE.optimizationLagrangeTermName), makeUntypedCrefIdent("dummyVarC"), "C"))%>; return 0; - #endif + >> + << + *res = &<%cref(makeUntypedCrefIdent(BackendDAE.optimizationLagrangeTermName))%>; + <%setLagrangeIndex%> >> let setInput = match simCode diff --git a/Compiler/Template/CodegenCFunctions.tpl b/Compiler/Template/CodegenCFunctions.tpl index c2633149374..1a45ff82b7f 100644 --- a/Compiler/Template/CodegenCFunctions.tpl +++ b/Compiler/Template/CodegenCFunctions.tpl @@ -4045,9 +4045,45 @@ template contextCref(ComponentRef cr, Context context, Text &auxFunction) >> else "_" + System.unquoteIdentifier(crefStr(cr)) ) + case JACOBIAN_CONTEXT(jacHT=SOME(_)) then jacCrefs(cr, context, 0) else cref(cr) end contextCref; +template contextCrefOld(ComponentRef cr, Context context, Text &auxFunction, Integer ix) + "Generates code for a component reference depending on which context we're in." +::= + match context + case FUNCTION_CONTEXT(__) + case PARALLEL_FUNCTION_CONTEXT(__) then + (match cr + case CREF_QUAL(identType = T_ARRAY(ty = T_COMPLEX(complexClassType = record_state))) then + let &preExp = buffer "" + let &varDecls = buffer "" + let rec_name = '<%underscorePath(ClassInf.getStateName(record_state))%>' + let recPtr = tempDecl(rec_name + "*", &varDecls) + let dimsLenStr = listLength(crefSubs(cr)) + let dimsValuesStr = (crefSubs(cr) |> INDEX(__) => daeDimensionExp(exp, context, &preExp, &varDecls, &auxFunction) ; separator=", ") + << + ((<%rec_name%>*)(generic_array_element_addr(&_<%ident%>, sizeof(<%rec_name%>), <%dimsLenStr%>, <%dimsValuesStr%>)))-><%contextCref(componentRef, context, &auxFunction)%> + >> + else "_" + System.unquoteIdentifier(crefStr(cr)) + ) + case JACOBIAN_CONTEXT(jacHT=SOME(_)) then jacCrefs(cr, context, ix) + else crefOld(cr, ix) +end contextCrefOld; + +template jacCrefs(ComponentRef cr, Context context, Integer ix) + "Generates code for jacobian variables." +::= + match context + case JACOBIAN_CONTEXT(jacHT=SOME(jacHT)) then + match simVarFromHT(cr, jacHT) + case SIMVAR(varKind=BackendDAE.JAC_VAR()) then 'jacobian->resultVars[<%index%>] /* <%escapeCComments(crefStrNoUnderscore(name))%> <%variabilityString(varKind)%> */' + case SIMVAR(varKind=BackendDAE.JAC_DIFF_VAR()) then 'jacobian->tmpVars[<%index%>] /* <%escapeCComments(crefStrNoUnderscore(name))%> <%variabilityString(varKind)%> */' + case SIMVAR(varKind=BackendDAE.SEED_VAR()) then 'jacobian->seedVars[<%index%>] /* <%escapeCComments(crefStrNoUnderscore(name))%> <%variabilityString(varKind)%> */' + case SIMVAR(index=-2) then crefOld(cr, ix) +end jacCrefs; + template contextCrefIsPre(ComponentRef cr, Context context, Text &auxFunction, Boolean isPre) "Generates code for a component reference depending on which context we're in." ::= @@ -4084,7 +4120,9 @@ end cref; ::= match cr case CREF_IDENT(ident = "xloc") then crefStr(cr) - case CREF_IDENT(ident = "time") then "data->localData[<%ix%>]->timeValue" + case CREF_IDENT(ident = "time") then 'data->localData[<%ix%>]->timeValue' + case CREF_IDENT(ident = "__OMC_DT") then "data->simulationInfo->inlineData->dt" + case CREF_IDENT(ident = "__HOM_LAMBDA") then "data->simulationInfo->lambda" case WILD(__) then '' else crefToCStr(cr, ix, false, false) end crefOld; @@ -4125,9 +4163,9 @@ template crefToCStr(ComponentRef cr, Integer ix, Boolean isPre, Boolean isStart) case SIMVAR(aliasvar=ALIAS(varName=varName)) then crefToCStr(varName, ix, isPre, isStart) case SIMVAR(aliasvar=NEGATEDALIAS(varName=varName), type_=T_BOOL()) then '!(<%crefToCStr(varName, ix, isPre, isStart)%>)' case SIMVAR(aliasvar=NEGATEDALIAS(varName=varName)) then '-(<%crefToCStr(varName, ix, isPre, isStart)%>)' - case SIMVAR(varKind=JAC_VAR()) - case SIMVAR(varKind=JAC_DIFF_VAR()) - case SIMVAR(varKind=SEED_VAR()) + case SIMVAR(varKind=JAC_VAR()) then 'parentJacobian->resultVars[<%index%>] /* <%escapeCComments(crefStrNoUnderscore(name))%> <%variabilityString(varKind)%> */' + case SIMVAR(varKind=JAC_DIFF_VAR()) then 'parentJacobian->tmpVars[<%index%>] /* <%escapeCComments(crefStrNoUnderscore(name))%> <%variabilityString(varKind)%> */' + case SIMVAR(varKind=SEED_VAR()) then 'parentJacobian->seedVars[<%index%>] /* <%escapeCComments(crefStrNoUnderscore(name))%> <%variabilityString(varKind)%> */' case SIMVAR(varKind=DAE_RESIDUAL_VAR()) case SIMVAR(varKind=DAE_AUX_VAR()) case SIMVAR(index=-2) @@ -5229,6 +5267,7 @@ template daeExpRelationSim(Exp exp, Context context, Text &preExp, match exp case rel as RELATION(__) then match context + case JACOBIAN_CONTEXT(__) case DAE_MODE_CONTEXT(__) case SIMULATION_CONTEXT(__) then match rel.optionExpisASUB diff --git a/Compiler/Template/SimCodeTV.mo b/Compiler/Template/SimCodeTV.mo index c783d81f300..f77725bdaef 100644 --- a/Compiler/Template/SimCodeTV.mo +++ b/Compiler/Template/SimCodeTV.mo @@ -280,6 +280,41 @@ package SimCodeVar end Causality; end SimCodeVar; +package HashTableCrefSimVar + + type Key = DAE.ComponentRef; + type Value = SimCodeVar.SimVar; + + type HashTableCrefFunctionsType = tuple; + type HashTable = tuple< + array>>, + tuple>>>, + Integer, + HashTableCrefFunctionsType + >; + + function FuncHashCref + input Key cr; + input Integer mod; + output Integer res; + end FuncHashCref; + + function FuncCrefEqual + input Key cr1; + input Key cr2; + output Boolean res; + end FuncCrefEqual; + + function FuncCrefStr + input Key cr; + output String res; + end FuncCrefStr; + + function FuncExpStr + input Value exp; + output String res; + end FuncExpStr; +end HashTableCrefSimVar; package SimCode @@ -308,6 +343,7 @@ package SimCode Integer maxColorCols; Integer jacobianIndex; Integer partitionIndex; + Option crefsHT; end JAC_MATRIX; end JacobianMatrix; @@ -532,6 +568,7 @@ package SimCode list sources; Integer indexLinearSystem; Integer nUnknowns; + Boolean partOfJac "if TRUE then this system is part of a jacobian matrix"; end LINEARSYSTEM; end LinearSystem; @@ -819,6 +856,7 @@ package SimCodeFunction record FUNCTION_CONTEXT end FUNCTION_CONTEXT; record JACOBIAN_CONTEXT + Option jacHT; end JACOBIAN_CONTEXT; record ALGLOOP_CONTEXT Boolean genInitialisation; @@ -1028,6 +1066,17 @@ package SimCodeUtil output SimCodeVar.SimVar outSimVar; end cref2simvar; + function simVarFromHT + input DAE.ComponentRef inCref; + input HashTableCrefSimVar.HashTable crefToSimVarHT; + output SimCodeVar.SimVar outSimVar; + end simVarFromHT; + + function createJacContext + input Option jacHT; + output SimCodeFunction.Context outContext; + end createJacContext; + function isModelTooBigForCSharpInOneFile input SimCode.SimCode simCode; output Boolean outIsTooBig; @@ -3331,6 +3380,13 @@ package ComponentReference input DAE.ComponentRef inCR; output DAE.ComponentRef outCR; end crefRemovePrePrefix; + + function createDifferentiatedCrefName + input DAE.ComponentRef inCref; + input DAE.ComponentRef inX; + input String inMatrixName; + output DAE.ComponentRef outCref; + end createDifferentiatedCrefName; end ComponentReference; package Expression diff --git a/Compiler/Util/Error.mo b/Compiler/Util/Error.mo index 92a8f01e6ba..51d43fae176 100644 --- a/Compiler/Util/Error.mo +++ b/Compiler/Util/Error.mo @@ -1144,6 +1144,9 @@ public constant Message SAVE_ENCRYPTED_CLASS_ERROR = MESSAGE(7022, SCRIPTING(), Util.gettext("Cannot save the encrypted class. Encrypted classes are read-only.")); public constant Message ACCESS_ENCRYPTED_PROTECTED_CONTENTS = MESSAGE(7023, SCRIPTING(), ERROR(), Util.gettext("Cannot access encrypted and protected class contents.")); +public constant Message INVALID_NONLINEAR_JACOBIAN_COMPONENT = MESSAGE(7024, TRANSLATION(), ERROR(), + Util.gettext("Jacobian %s contains non-linear components. This indicates a singular system or internal generation errors.")); + protected diff --git a/SimulationRuntime/c/linearization/linearize.cpp b/SimulationRuntime/c/linearization/linearize.cpp index 39657dd1b92..8f21674f973 100644 --- a/SimulationRuntime/c/linearization/linearize.cpp +++ b/SimulationRuntime/c/linearization/linearize.cpp @@ -288,37 +288,38 @@ int functionJacBD_num(DATA* data, threadData_t *threadData, double *matrixB, dou int functionJacA(DATA* data, threadData_t *threadData, double* jac){ const int index = data->callback->INDEX_JAC_A; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); unsigned int i,j,k; k = 0; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeCols; i++) + for(i=0; i < jacobian->sizeCols; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1.0; + jacobian->seedVars[i] = 1.0; if(ACTIVE_STREAM(LOG_JAC)) { printf("Caluculate one col:\n"); - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) + for(j=0; j < jacobian->sizeCols;j++) { - infoStreamPrint(LOG_JAC,0,"seed: data->simulationInfo->analyticJacobians[index].seedVars[%d]= %f",j,data->simulationInfo->analyticJacobians[index].seedVars[j]); + infoStreamPrint(LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]); } } - data->callback->functionJacA_column(data, threadData); + data->callback->functionJacA_column(data, threadData, jacobian, NULL); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeRows; j++) + for(j = 0; j < jacobian->sizeRows; j++) { - jac[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j]; - infoStreamPrint(LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,data->simulationInfo->analyticJacobians[index].resultVars[j]); + jac[k++] = jacobian->resultVars[j]; + infoStreamPrint(LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]); } - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0.0; + jacobian->seedVars[i] = 0.0; } if(ACTIVE_STREAM(LOG_JAC)) { infoStreamPrint(LOG_JAC,0,"Print jac:"); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows;i++) + for(i=0; i < jacobian->sizeRows;i++) { - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) { - printf("% .5e ",jac[i+j*data->simulationInfo->analyticJacobians[index].sizeCols]); + for(j=0; j < jacobian->sizeCols;j++) { + printf("% .5e ",jac[i+j*jacobian->sizeCols]); } printf("\n"); } @@ -329,37 +330,39 @@ int functionJacA(DATA* data, threadData_t *threadData, double* jac){ int functionJacB(DATA* data, threadData_t *threadData, double* jac){ const int index = data->callback->INDEX_JAC_B; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); + unsigned int i,j,k; k = 0; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeCols; i++) + for(i=0; i < jacobian->sizeCols; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1.0; + jacobian->seedVars[i] = 1.0; if(ACTIVE_STREAM(LOG_JAC)) { printf("Caluculate one col:\n"); - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) + for(j=0; j < jacobian->sizeCols;j++) { - infoStreamPrint(LOG_JAC,0,"seed: data->simulationInfo->analyticJacobians[index].seedVars[%d]= %f",j,data->simulationInfo->analyticJacobians[index].seedVars[j]); + infoStreamPrint(LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]); } } - data->callback->functionJacB_column(data, threadData); + data->callback->functionJacB_column(data, threadData, jacobian, NULL); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeRows; j++) + for(j = 0; j < jacobian->sizeRows; j++) { - jac[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j]; - infoStreamPrint(LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,data->simulationInfo->analyticJacobians[index].resultVars[j]); + jac[k++] = jacobian->resultVars[j]; + infoStreamPrint(LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]); } - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0.0; + jacobian->seedVars[i] = 0.0; } if(ACTIVE_STREAM(LOG_JAC)) { infoStreamPrint(LOG_JAC, 0, "Print jac:"); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows;i++) + for(i=0; i < jacobian->sizeRows;i++) { - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) - printf("% .5e ",jac[i+j*data->simulationInfo->analyticJacobians[index].sizeCols]); + for(j=0; j < jacobian->sizeCols;j++) + printf("% .5e ",jac[i+j*jacobian->sizeCols]); printf("\n"); } } @@ -369,35 +372,36 @@ int functionJacB(DATA* data, threadData_t *threadData, double* jac){ int functionJacC(DATA* data, threadData_t *threadData, double* jac){ const int index = data->callback->INDEX_JAC_C; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); unsigned int i,j,k; k = 0; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeCols; i++) + for(i=0; i < jacobian->sizeCols; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1.0; + jacobian->seedVars[i] = 1.0; if(ACTIVE_STREAM(LOG_JAC)) { - printf("Calculate one col:\n"); - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) - infoStreamPrint(LOG_JAC,0,"seed: data->simulationInfo->analyticJacobians[index].seedVars[%d]= %f",j,data->simulationInfo->analyticJacobians[index].seedVars[j]); + printf("Caluculate one col:\n"); + for(j=0; j < jacobian->sizeCols;j++) + infoStreamPrint(LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]); } - data->callback->functionJacC_column(data, threadData); + data->callback->functionJacC_column(data, threadData, jacobian, NULL); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeRows; j++) + for(j = 0; j < jacobian->sizeRows; j++) { - jac[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j]; - infoStreamPrint(LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,data->simulationInfo->analyticJacobians[index].resultVars[j]); + jac[k++] = jacobian->resultVars[j]; + infoStreamPrint(LOG_JAC,0,"write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]); } - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0.0; + jacobian->seedVars[i] = 0.0; } if(ACTIVE_STREAM(LOG_JAC)) { infoStreamPrint(LOG_JAC, 0, "Print jac:"); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows;i++) + for(i=0; i < jacobian->sizeRows;i++) { - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) - printf("% .5e ",jac[i+j*data->simulationInfo->analyticJacobians[index].sizeCols]); + for(j=0; j < jacobian->sizeCols;j++) + printf("% .5e ",jac[i+j*jacobian->sizeCols]); printf("\n"); } } @@ -407,36 +411,37 @@ int functionJacC(DATA* data, threadData_t *threadData, double* jac){ int functionJacD(DATA* data, threadData_t *threadData, double* jac){ const int index = data->callback->INDEX_JAC_D; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); unsigned int i,j,k; k = 0; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeCols; i++) + for(i=0; i < jacobian->sizeCols; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1.0; + jacobian->seedVars[i] = 1.0; if(ACTIVE_STREAM(LOG_JAC)) { - printf("Calculate one col:\n"); - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) { - infoStreamPrint(LOG_JAC,0,"seed: data->simulationInfo->analyticJacobians[index].seedVars[%d]= %f",j,data->simulationInfo->analyticJacobians[index].seedVars[j]); + printf("Caluculate one col:\n"); + for(j=0; j < jacobian->sizeCols;j++) { + infoStreamPrint(LOG_JAC,0,"seed: jacobian->seedVars[%d]= %f",j,jacobian->seedVars[j]); } } - data->callback->functionJacD_column(data, threadData); + data->callback->functionJacD_column(data, threadData, jacobian, NULL); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeRows; j++) + for(j = 0; j < jacobian->sizeRows; j++) { - jac[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j]; - infoStreamPrint(LOG_JAC,0, "write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,data->simulationInfo->analyticJacobians[index].resultVars[j]); + jac[k++] = jacobian->resultVars[j]; + infoStreamPrint(LOG_JAC,0, "write in jac[%d]-[%d,%d]=%g from row[%d]=%g",k-1,i,j,jac[k-1],i,jacobian->resultVars[j]); } - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0.0; + jacobian->seedVars[i] = 0.0; } if(ACTIVE_STREAM(LOG_JAC)) { infoStreamPrint(LOG_JAC, 0, "Print jac:"); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows;i++) + for(i=0; i < jacobian->sizeRows;i++) { - for(j=0; j < data->simulationInfo->analyticJacobians[index].sizeCols;j++) - printf("% .5e ",jac[i+j*data->simulationInfo->analyticJacobians[index].sizeCols]); + for(j=0; j < jacobian->sizeCols;j++) + printf("% .5e ",jac[i+j*jacobian->sizeCols]); printf("\n"); } } @@ -508,22 +513,26 @@ int linearize(DATA* data, threadData_t *threadData) if (data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sizeTmpVars > 0){ /* Retrieve symbolic Jacobian */ /* Determine Matrix A */ - if(!data->callback->initialAnalyticJacobianA(data, threadData)){ + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]); + if(!data->callback->initialAnalyticJacobianA(data, threadData, jacobian)){ assertStreamPrint(threadData,0==functionJacA(data, threadData, matrixA),"Error, can not get Matrix A "); } /* Determine Matrix B */ - if(!data->callback->initialAnalyticJacobianB(data, threadData)){ + jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_B]); + if(!data->callback->initialAnalyticJacobianB(data, threadData, jacobian)){ assertStreamPrint(threadData,0==functionJacB(data, threadData, matrixB),"Error, can not get Matrix B "); } /* Determine Matrix C */ - if(!data->callback->initialAnalyticJacobianC(data, threadData)){ + jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_C]); + if(!data->callback->initialAnalyticJacobianC(data, threadData, jacobian)){ assertStreamPrint(threadData,0==functionJacC(data, threadData, matrixC),"Error, can not get Matrix C "); } /* Determine Matrix D */ - if(!data->callback->initialAnalyticJacobianD(data, threadData)){ + jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_D]); + if(!data->callback->initialAnalyticJacobianD(data, threadData, jacobian)){ assertStreamPrint(threadData,0==functionJacD(data, threadData, matrixD),"Error, can not get Matrix D "); } } diff --git a/SimulationRuntime/c/openmodelica_func.h b/SimulationRuntime/c/openmodelica_func.h index ee671146a5d..13131a37c66 100644 --- a/SimulationRuntime/c/openmodelica_func.h +++ b/SimulationRuntime/c/openmodelica_func.h @@ -242,18 +242,18 @@ const int INDEX_JAC_D; * Return-value 0: jac is present * Return-value 1: jac is not present */ -int (*initialAnalyticJacobianA)(void* data, threadData_t *threadData); -int (*initialAnalyticJacobianB)(void* data, threadData_t *threadData); -int (*initialAnalyticJacobianC)(void* data, threadData_t *threadData); -int (*initialAnalyticJacobianD)(void* data, threadData_t *threadData); +int (*initialAnalyticJacobianA)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian); +int (*initialAnalyticJacobianB)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian); +int (*initialAnalyticJacobianC)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian); +int (*initialAnalyticJacobianD)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian); /* * These functions calculate specific jacobian column. */ -int (*functionJacA_column)(void* data, threadData_t *threadData); -int (*functionJacB_column)(void* data, threadData_t *threadData); -int (*functionJacC_column)(void* data, threadData_t *threadData); -int (*functionJacD_column)(void* data, threadData_t *threadData); +int (*functionJacA_column)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian, ANALYTIC_JACOBIAN* parentJacobian); +int (*functionJacB_column)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian, ANALYTIC_JACOBIAN* parentJacobian); +int (*functionJacC_column)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian, ANALYTIC_JACOBIAN* parentJacobian); +int (*functionJacD_column)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian, ANALYTIC_JACOBIAN* parentJacobian); /*#endif*/ @@ -332,8 +332,8 @@ void (*read_input_fmu)(MODEL_DATA* modelData, SIMULATION_INFO* simulationData); /* * FMU continuous partial derivative functions */ -int (*initialPartialFMIDER)(void* data, threadData_t *threadData); -int (*functionJacFMIDER_column)(void* data, threadData_t *threadData); +int (*initialPartialFMIDER)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian); +int (*functionJacFMIDER_column)(void* data, threadData_t *threadData, ANALYTIC_JACOBIAN* thisJacobian, ANALYTIC_JACOBIAN* parentJacobian); const int INDEX_JAC_FMIDER; }; diff --git a/SimulationRuntime/c/optimization/DataManagement/DerStructure.c b/SimulationRuntime/c/optimization/DataManagement/DerStructure.c index 32e9d17d0c3..825c458a45b 100644 --- a/SimulationRuntime/c/optimization/DataManagement/DerStructure.c +++ b/SimulationRuntime/c/optimization/DataManagement/DerStructure.c @@ -74,10 +74,10 @@ inline void allocate_der_struct(OptDataStructure *s, OptDataDim * dim, DATA* dat s->indexABCD[3] = data->callback->INDEX_JAC_C; s->indexABCD[4] = data->callback->INDEX_JAC_D; s->matrix[0] = (modelica_boolean)0; - s->matrix[1] = (modelica_boolean)(data->callback->initialAnalyticJacobianA((void*) data, threadData) == 0); - s->matrix[2] = (modelica_boolean)(data->callback->initialAnalyticJacobianB((void*) data, threadData) == 0); - s->matrix[3] = (modelica_boolean)(data->callback->initialAnalyticJacobianC((void*) data, threadData) == 0); - s->matrix[4] = (modelica_boolean)(data->callback->initialAnalyticJacobianD((void*) data, threadData) == 0); + s->matrix[1] = (modelica_boolean)(data->callback->initialAnalyticJacobianA((void*) data, threadData, &(data->simulationInfo->analyticJacobians[s->indexABCD[1]])) == 0); + s->matrix[2] = (modelica_boolean)(data->callback->initialAnalyticJacobianB((void*) data, threadData, &(data->simulationInfo->analyticJacobians[s->indexABCD[2]])) == 0); + s->matrix[3] = (modelica_boolean)(data->callback->initialAnalyticJacobianC((void*) data, threadData, &(data->simulationInfo->analyticJacobians[s->indexABCD[3]])) == 0); + s->matrix[4] = (modelica_boolean)(data->callback->initialAnalyticJacobianD((void*) data, threadData, &(data->simulationInfo->analyticJacobians[s->indexABCD[4]])) == 0); dim->nJderx = 0; dim->nJfderx = 0; diff --git a/SimulationRuntime/c/optimization/DataManagement/MoveData.c b/SimulationRuntime/c/optimization/DataManagement/MoveData.c index 7b4778e9a90..766e2ed2818 100644 --- a/SimulationRuntime/c/optimization/DataManagement/MoveData.c +++ b/SimulationRuntime/c/optimization/DataManagement/MoveData.c @@ -817,15 +817,16 @@ void diffSynColoredOptimizerSystem(OptData *optData, modelica_real **J, const in int i,j,l,ii, ll; const int h_index = optData->s.indexABCD[index]; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]); const long double * scaldt = optData->bounds.scaldt[m]; - const unsigned int * const cC = data->simulationInfo->analyticJacobians[h_index].sparsePattern.colorCols; - const unsigned int * const lindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.leadindex; - const int nx = data->simulationInfo->analyticJacobians[h_index].sizeCols; - const int Cmax = data->simulationInfo->analyticJacobians[h_index].sparsePattern.maxColors + 1; + const unsigned int * const cC = jacobian->sparsePattern.colorCols; + const unsigned int * const lindex = jacobian->sparsePattern.leadindex; + const int nx = jacobian->sizeCols; + const int Cmax = jacobian->sparsePattern.maxColors + 1; const int dnx = optData->dim.nx; const int dnxnc = optData->dim.nJ; - const modelica_real * const resultVars = data->simulationInfo->analyticJacobians[h_index].resultVars; - const unsigned int * const sPindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.index; + const modelica_real * const resultVars = jacobian->resultVars; + const unsigned int * const sPindex = jacobian->sparsePattern.index; long double scalb = optData->bounds.scalb[m][n]; const int * index_J = (index == 3)? optData->s.indexJ3 : optData->s.indexJ2; @@ -837,12 +838,12 @@ void diffSynColoredOptimizerSystem(OptData *optData, modelica_real **J, const in setContext(data, &(data->localData[0]->timeValue), CONTEXT_SYM_JACOBIAN); for(i = 1; i < Cmax; ++i){ - data->simulationInfo->analyticJacobians[h_index].seedVars = sV[i]; + jacobian->seedVars = sV[i]; if(index == 2){ - data->callback->functionJacB_column(data, threadData); + data->callback->functionJacB_column(data, threadData, jacobian, NULL); }else if(index == 3){ - data->callback->functionJacC_column(data, threadData); + data->callback->functionJacC_column(data, threadData, jacobian, NULL); }else assert(0); @@ -878,12 +879,13 @@ void diffSynColoredOptimizerSystemF(OptData *optData, modelica_real **J){ int i,j,l,ii, ll; const int index = 4; const int h_index = optData->s.indexABCD[index]; - const unsigned int * const cC = data->simulationInfo->analyticJacobians[h_index].sparsePattern.colorCols; - const unsigned int * const lindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.leadindex; - const int nx = data->simulationInfo->analyticJacobians[h_index].sizeCols; - const int Cmax = data->simulationInfo->analyticJacobians[h_index].sparsePattern.maxColors + 1; - const modelica_real * const resultVars = data->simulationInfo->analyticJacobians[h_index].resultVars; - const unsigned int * const sPindex = data->simulationInfo->analyticJacobians[h_index].sparsePattern.index; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[h_index]); + const unsigned int * const cC = jacobian->sparsePattern.colorCols; + const unsigned int * const lindex = jacobian->sparsePattern.leadindex; + const int nx = jacobian->sizeCols; + const int Cmax = jacobian->sparsePattern.maxColors + 1; + const modelica_real * const resultVars = jacobian->resultVars; + const unsigned int * const sPindex = jacobian->sparsePattern.index; modelica_real **sV = optData->s.seedVec[index]; @@ -891,9 +893,9 @@ void diffSynColoredOptimizerSystemF(OptData *optData, modelica_real **J){ setContext(data, &(data->localData[0]->timeValue), CONTEXT_SYM_JACOBIAN); for(i = 1; i < Cmax; ++i){ - data->simulationInfo->analyticJacobians[h_index].seedVars = sV[i]; + jacobian->seedVars = sV[i]; - data->callback->functionJacD_column(data, threadData); + data->callback->functionJacD_column(data, threadData, jacobian, NULL); increaseJacContext(data); diff --git a/SimulationRuntime/c/simulation/solver/dassl.c b/SimulationRuntime/c/simulation/solver/dassl.c index a59bb252bf7..48431cdea98 100644 --- a/SimulationRuntime/c/simulation/solver/dassl.c +++ b/SimulationRuntime/c/simulation/solver/dassl.c @@ -306,7 +306,8 @@ int dassl_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo, dasslData->dasslJacobian == COLOREDSYMJAC || dasslData->dasslJacobian == SYMJAC) { - if (data->callback->initialAnalyticJacobianA(data, threadData)) + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]); + if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian)) { infoStreamPrint(LOG_STDOUT, 0, "Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal."); dasslData->dasslJacobian = INTERNALNUMJAC; @@ -854,38 +855,40 @@ int jacA_symColored(double *t, double *y, double *yprime, double *delta, double threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2]; const int index = data->callback->INDEX_JAC_A; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); + unsigned int i,j,l,k,ii; /* set symbolical jacobian to reuse the matrix A and the factorization * in the Linear loops of functionJacA_column */ setContext(data, t, CONTEXT_SYM_JACOBIAN); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; - data->callback->functionJacA_column(data, threadData); + data->callback->functionJacA_column(data, threadData, jacobian, NULL); increaseJacContext(data); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l; - matrixA[k] = data->simulationInfo->analyticJacobians[index].resultVars[l]; + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; + matrixA[k] = jacobian->resultVars[l]; ii++; }; } } - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) data->simulationInfo->analyticJacobians[index].seedVars[ii] = 0; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) jacobian->seedVars[ii] = 0; } @@ -905,7 +908,10 @@ int jacA_sym(double *t, double *y, double *yprime, double *delta, double *matrix DATA* data = (DATA*)(void*)((double**)rpar)[0]; DASSL_DATA* dasslData = (DASSL_DATA*)(void*)((double**)rpar)[1]; threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2]; + const int index = data->callback->INDEX_JAC_A; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); + unsigned int i,j,k; /* set symbolical jacobian to reuse the matrix A and the factorization @@ -913,20 +919,20 @@ int jacA_sym(double *t, double *y, double *yprime, double *delta, double *matrix setContext(data, t, CONTEXT_SYM_JACOBIAN); k = 0; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeCols; i++) + for(i=0; i < jacobian->sizeCols; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1.0; + jacobian->seedVars[i] = 1.0; - data->callback->functionJacA_column(data, threadData); + data->callback->functionJacA_column(data, threadData, jacobian, NULL); increaseJacContext(data); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeRows; j++) + for(j = 0; j < jacobian->sizeRows; j++) { - matrixA[k++] = data->simulationInfo->analyticJacobians[index].resultVars[j]; + matrixA[i*jacobian->sizeCols+j] = jacobian->resultVars[j]; } - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0.0; + jacobian->seedVars[i] = 0.0; } TRACE_POP @@ -1004,6 +1010,8 @@ int jacA_numColored(double *t, double *y, double *yprime, double *delta, double threadData_t *threadData = (threadData_t*)(void*)((double**)rpar)[2]; const int index = data->callback->INDEX_JAC_A; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); + double delta_h = numericalDifferentiationDeltaXsolver; double delta_hhh; int ires; @@ -1016,11 +1024,11 @@ int jacA_numColored(double *t, double *y, double *yprime, double *delta, double /* set context for the start values extrapolation of non-linear algebraic loops */ setContext(data, t, CONTEXT_JACOBIAN); - for(i = 0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i = 0; i < jacobian->sparsePattern.maxColors; i++) { - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) + for(ii=0; ii < jacobian->sizeCols; ii++) { - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) { delta_hhh = *h * yprime[ii]; delta_hh[ii] = delta_h * fmax(fmax(fabs(y[ii]),fabs(delta_hhh)),fabs(1./wt[ii])); @@ -1038,15 +1046,15 @@ int jacA_numColored(double *t, double *y, double *yprime, double *delta, double increaseJacContext(data); - for(ii = 0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) + for(ii = 0; ii < jacobian->sizeCols; ii++) { - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) { - j = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii]; - while(j < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii+1]) + j = jacobian->sparsePattern.leadindex[ii]; + while(j < jacobian->sparsePattern.leadindex[ii+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[j]; - k = l + ii*data->simulationInfo->analyticJacobians[index].sizeRows; + l = jacobian->sparsePattern.index[j]; + k = l + ii*jacobian->sizeRows; matrixA[k] = (dasslData->newdelta[l] - delta[l]) * delta_hh[ii]; j++; }; diff --git a/SimulationRuntime/c/simulation/solver/ida_solver.c b/SimulationRuntime/c/simulation/solver/ida_solver.c index c544629cf26..6df03401de1 100644 --- a/SimulationRuntime/c/simulation/solver/ida_solver.c +++ b/SimulationRuntime/c/simulation/solver/ida_solver.c @@ -440,7 +440,8 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo idaData->jacobianMethod == COLOREDSYMJAC || idaData->jacobianMethod == SYMJAC)) { - if (data->callback->initialAnalyticJacobianA(data, threadData)) + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A]); + if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian)) { infoStreamPrint(LOG_STDOUT, 0, "Jacobian or SparsePattern is not generated or failed to initialize! Switch back to normal."); idaData->jacobianMethod = INTERNALNUMJAC; @@ -1444,7 +1445,7 @@ int jacColoredSymbolicalDense(double tt, N_Vector yy, N_Vector yp, N_Vector rr, } } - data->callback->functionJacA_column(data, threadData); + data->callback->functionJacA_column(data, threadData, jacData, NULL); increaseJacContext(data); for(ii = 0; ii < idaData->N; ii++) @@ -1724,7 +1725,7 @@ jacColoredSymbolicalSparse(double tt, N_Vector yy, N_Vector yp, N_Vector rr, Sls } } - data->callback->functionJacA_column(data, threadData); + data->callback->functionJacA_column(data, threadData, jacData, NULL); increaseJacContext(data); for(ii = 0; ii < idaData->N; ii++) diff --git a/SimulationRuntime/c/simulation/solver/kinsolSolver.c b/SimulationRuntime/c/simulation/solver/kinsolSolver.c index 541df69a317..a27fd237662 100644 --- a/SimulationRuntime/c/simulation/solver/kinsolSolver.c +++ b/SimulationRuntime/c/simulation/solver/kinsolSolver.c @@ -554,7 +554,7 @@ int nlsSparseSymJac(N_Vector vecX, N_Vector vecFX, SlsMat Jac, void *userData, N analyticJacobian->seedVars[ii] = 1.0; } } - ((nlsData->analyticalJacobianColumn))(data, threadData); + ((nlsData->analyticalJacobianColumn))(data, threadData, analyticJacobian, NULL); for(ii = 0; ii < kinsolData->size; ii++) { diff --git a/SimulationRuntime/c/simulation/solver/linearSolverKlu.c b/SimulationRuntime/c/simulation/solver/linearSolverKlu.c index 1086808427e..d7b704a689b 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverKlu.c +++ b/SimulationRuntime/c/simulation/solver/linearSolverKlu.c @@ -123,24 +123,27 @@ int getAnalyticalJacobian(DATA* data, threadData_t *threadData, int sysNumber) LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); + ANALYTIC_JACOBIAN* parentJacobian = systemData->parentJacobian; + int nth = 0; - int nnz = data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros; + int nnz = jacobian->sparsePattern.numberOfNoneZeros; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows; i++) + for(i=0; i < jacobian->sizeRows; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1; + jacobian->seedVars[i] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, parentJacobian); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - systemData->setAElement(i, l, -data->simulationInfo->analyticJacobians[index].resultVars[l], nth, (void*) systemData, threadData); + l = jacobian->sparsePattern.index[ii]; + systemData->setAElement(i, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData); nth++; ii++; }; @@ -148,7 +151,7 @@ int getAnalyticalJacobian(DATA* data, threadData_t *threadData, int sysNumber) }; /* de-activate seed variable for the corresponding color */ - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0; + jacobian->seedVars[i] = 0; } return 0; @@ -174,7 +177,7 @@ static int residual_wrapper(double* x, double* f, void** data, int sysNumber) * author: wbraun */ int -solveKlu(DATA *data, threadData_t *threadData, int sysNumber) +solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); @@ -214,7 +217,7 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber) } /* calculate vector b (rhs) */ - memcpy(solverData->work, systemData->x, sizeof(double)*solverData->n_row); + memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row); residual_wrapper(solverData->work, systemData->b, dataAndThreadData, sysNumber); } tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock)); @@ -225,7 +228,7 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber) { infoStreamPrint(LOG_LS_V, 1, "Old solution x:"); for(i = 0; i < solverData->n_row; ++i) - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); messageClose(LOG_LS_V); infoStreamPrint(LOG_LS_V, 1, "Matrix A n_rows = %d", solverData->n_row); @@ -292,13 +295,13 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber) if (1 == systemData->method){ /* take the solution */ for(i = 0; i < solverData->n_row; ++i) - systemData->x[i] += systemData->b[i]; + aux_x[i] += systemData->b[i]; /* update inner equations */ - residual_wrapper(systemData->x, solverData->work, dataAndThreadData, sysNumber); + residual_wrapper(aux_x, solverData->work, dataAndThreadData, sysNumber); } else { /* the solution is automatically in x */ - memcpy(systemData->x, systemData->b, sizeof(double)*systemData->size); + memcpy(aux_x, systemData->b, sizeof(double)*systemData->size); } if (ACTIVE_STREAM(LOG_LS_V)) @@ -307,7 +310,7 @@ solveKlu(DATA *data, threadData_t *threadData, int sysNumber) infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar); for(i = 0; i < systemData->size; ++i) - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); messageClose(LOG_LS_V); } diff --git a/SimulationRuntime/c/simulation/solver/linearSolverKlu.h b/SimulationRuntime/c/simulation/solver/linearSolverKlu.h index e6c8f7859ab..373e15c56ab 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverKlu.h +++ b/SimulationRuntime/c/simulation/solver/linearSolverKlu.h @@ -62,7 +62,7 @@ typedef struct DATA_KLU int allocateKluData(int n_row, int n_col, int nz, void **data); int freeKluData(void **data); -int solveKlu(DATA *data, threadData_t *threadData, int sysNumber); +int solveKlu(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); #endif #endif diff --git a/SimulationRuntime/c/simulation/solver/linearSolverLapack.c b/SimulationRuntime/c/simulation/solver/linearSolverLapack.c index adb731475e5..83dffc6f3a9 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverLapack.c +++ b/SimulationRuntime/c/simulation/solver/linearSolverLapack.c @@ -108,34 +108,35 @@ int getAnalyticalJacobianLapack(DATA* data, threadData_t *threadData, double* ja LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[currentSys]); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); memset(jac, 0, (systemData->size)*(systemData->size)*sizeof(double)); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { /* activate seed variable for the corresponding color */ - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, systemData->parentJacobian); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l; - jac[k] = -data->simulationInfo->analyticJacobians[index].resultVars[l]; + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; + jac[k] = -jacobian->resultVars[l]; ii++; }; } /* de-activate seed variable for the corresponding color */ - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[j]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[j] = 0; + if(jacobian->sparsePattern.colorCols[j]-1 == i) + jacobian->seedVars[j] = 0; } } @@ -160,7 +161,7 @@ static int wrapper_fvec_lapack(_omc_vector* x, _omc_vector* f, int* iflag, void* * * \author wbraun */ -int solveLapack(DATA *data, threadData_t *threadData, int sysNumber) +int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; int i, iflag = 1; @@ -183,7 +184,7 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber) /* set data */ - _omc_setVectorData(solverData->x, systemData->x); + _omc_setVectorData(solverData->x, aux_x); _omc_setVectorData(solverData->b, systemData->b); _omc_setMatrixData(solverData->A, systemData->A); @@ -307,7 +308,7 @@ int solveLapack(DATA *data, threadData_t *threadData, int sysNumber) infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar); for(i = 0; i < systemData->size; ++i) { - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %.15g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %.15g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); } messageClose(LOG_LS_V); diff --git a/SimulationRuntime/c/simulation/solver/linearSolverLapack.h b/SimulationRuntime/c/simulation/solver/linearSolverLapack.h index 3446af34777..b1ab4fc5306 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverLapack.h +++ b/SimulationRuntime/c/simulation/solver/linearSolverLapack.h @@ -54,7 +54,7 @@ typedef struct DATA_LAPACK int allocateLapackData(int size, void **data); int freeLapackData(void **data); -int solveLapack(DATA *data, threadData_t *threadData, int sysNumber); +int solveLapack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); #endif diff --git a/SimulationRuntime/c/simulation/solver/linearSolverLis.c b/SimulationRuntime/c/simulation/solver/linearSolverLis.c index b1475054596..e67172005c8 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverLis.c +++ b/SimulationRuntime/c/simulation/solver/linearSolverLis.c @@ -143,31 +143,33 @@ int getAnalyticalJacobianLis(DATA* data, threadData_t *threadData, int sysNumber LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); + int nth = 0; - int nnz = data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros; + int nnz = jacobian->sparsePattern.numberOfNoneZeros; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows; i++) + for(i=0; i < jacobian->sizeRows; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1; + jacobian->seedVars[i] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, systemData->parentJacobian); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j-1]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j-1]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - /*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -data->simulationInfo->analyticJacobians[index].resultVars[l]); */ - systemData->setAElement(i, l, -data->simulationInfo->analyticJacobians[index].resultVars[l], nth, (void*) systemData, threadData); + l = jacobian->sparsePattern.index[ii]; + /*infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -jacobian->resultVars[l]); */ + systemData->setAElement(i, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData); nth++; ii++; }; } } - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0; + jacobian->seedVars[i] = 0; } return 0; @@ -193,7 +195,7 @@ static int wrapper_fvec_lis(double* x, double* f, void** data, int sysNumber) */ int -solveLis(DATA *data, threadData_t *threadData, int sysNumber) +solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); @@ -210,7 +212,7 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber) /* set old values as start value for the iteration */ for(i=0; ix[i], solverData->x); + err = lis_vector_set_value(LIS_INS_VALUE, i, aux_x[i], solverData->x); } rt_ext_tp_tick(&(solverData->timeClock)); @@ -237,7 +239,7 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber) lis_matrix_assemble(solverData->A); /* calculate vector b (rhs) */ - memcpy(solverData->work, systemData->x, sizeof(double)*solverData->n_row); + memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row); wrapper_fvec_lis(solverData->work, systemData->b, dataAndThreadData, sysNumber); /* set b vector */ for(i=0; imethod){ /* take the solution */ - lis_vector_get_values(solverData->x, 0, solverData->n_row, systemData->x); + lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x); for(i = 0; i < solverData->n_row; ++i) - systemData->x[i] += solverData->work[i]; + aux_x[i] += solverData->work[i]; /* update inner equations */ - wrapper_fvec_lis(systemData->x, solverData->work, dataAndThreadData, sysNumber); + wrapper_fvec_lis(aux_x, solverData->work, dataAndThreadData, sysNumber); } else { /* write solution */ - lis_vector_get_values(solverData->x, 0, solverData->n_row, systemData->x); + lis_vector_get_values(solverData->x, 0, solverData->n_row, aux_x); } if (ACTIVE_STREAM(LOG_LS_V)) @@ -300,7 +302,7 @@ solveLis(DATA *data, threadData_t *threadData, int sysNumber) infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar); for(i = 0; i < systemData->size; ++i) - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); messageClose(LOG_LS_V); } diff --git a/SimulationRuntime/c/simulation/solver/linearSolverLis.h b/SimulationRuntime/c/simulation/solver/linearSolverLis.h index 7e15faa7e3d..3d5849c0037 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverLis.h +++ b/SimulationRuntime/c/simulation/solver/linearSolverLis.h @@ -56,7 +56,7 @@ typedef struct DATA_LIS int allocateLisData(int n_row, int n_col, int nz, void **data); int freeLisData(void **data); -int solveLis(DATA *data, threadData_t *threadData, int sysNumber); +int solveLis(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); #endif diff --git a/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c b/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c index 419fcf100bb..e17f88b482a 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c +++ b/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.c @@ -324,33 +324,34 @@ int getAnalyticalJacobianTotalPivot(DATA* data, threadData_t *threadData, double LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[currentSys]); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); memset(jac, 0, (systemData->size)*(systemData->size)*sizeof(double)); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { /* activate seed variable for the corresponding color */ - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, systemData->parentJacobian); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l; - jac[k] = data->simulationInfo->analyticJacobians[index].resultVars[l]; + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; + jac[k] = jacobian->resultVars[l]; ii++; } } /* de-activate seed variable for the corresponding color */ - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[j]-1 == i) { - data->simulationInfo->analyticJacobians[index].seedVars[j] = 0; + if(jacobian->sparsePattern.colorCols[j]-1 == i) { + jacobian->seedVars[j] = 0; } } @@ -381,7 +382,7 @@ static int wrapper_fvec_totalpivot(double* x, double* f, void** data, int sysNum * * \author bbachmann */ -int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber) +int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; int i, j; @@ -404,9 +405,8 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber) eqSystemNumber, (int) systemData->size, data->localData[0]->timeValue); - debugVectorDoubleLS(LOG_LS_V,"SCALING",systemData->nominal,n); - debugVectorDoubleLS(LOG_LS_V,"Old VALUES",systemData->x,n); + debugVectorDoubleLS(LOG_LS_V,"Old VALUES",aux_x,n); rt_ext_tp_tick(&(solverData->timeClock)); if (0 == systemData->method){ @@ -431,7 +431,7 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber) assertStreamPrint(threadData, 1, "jacobian function pointer is invalid" ); } /* calculate vector b (rhs) -> -b is last column of matrix Ab */ - wrapper_fvec_totalpivot(systemData->x, solverData->Ab + n*n, dataAndThreadData, sysNumber); + wrapper_fvec_totalpivot(aux_x, solverData->Ab + n*n, dataAndThreadData, sysNumber); } tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock)); systemData->jacobianTime += tmpJacEvalTime; @@ -451,11 +451,11 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber) debugVectorDoubleLS(LOG_LS_V, "SOLUTION:", solverData->x, n+1); if (1 == systemData->method){ /* add the solution to old solution vector*/ - vecAddLS(n, systemData->x, solverData->x, systemData->x); - wrapper_fvec_totalpivot(systemData->x, solverData->b, dataAndThreadData, sysNumber); + vecAddLS(n, aux_x, solverData->x, aux_x); + wrapper_fvec_totalpivot(aux_x, solverData->b, dataAndThreadData, sysNumber); } else { /* take the solution */ - vecCopyLS(n, solverData->x, systemData->x); + vecCopyLS(n, solverData->x, aux_x); } if (ACTIVE_STREAM(LOG_LS_V)){ @@ -463,7 +463,7 @@ int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber) infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar); for(i=0; isize; ++i) { - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); } messageClose(LOG_LS_V); } diff --git a/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h b/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h index 3e2bdc4c837..38242031e57 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h +++ b/SimulationRuntime/c/simulation/solver/linearSolverTotalPivot.h @@ -53,7 +53,7 @@ typedef struct DATA_TOTALPIVOT int allocateTotalPivotData(int size, void** data); int freeTotalPivotData(void** data); -int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber); +int solveTotalPivot(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); #endif diff --git a/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c b/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c index aad25cb6b1e..a16edc9b2d6 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c +++ b/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.c @@ -50,7 +50,7 @@ void printMatrixCSC(int* Ap, int* Ai, double* Ax, int n); void printMatrixCSR(int* Ap, int* Ai, double* Ax, int n); -int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData); +int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x); /*! \fn allocate memory for linear system solver UmfPack * @@ -137,25 +137,27 @@ int getAnalyticalJacobianUmfPack(DATA* data, threadData_t *threadData, int sysNu LINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->linearSystemData[sysNumber]); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); + int nth = 0; - int nnz = data->simulationInfo->analyticJacobians[index].sparsePattern.numberOfNoneZeros; + int nnz = jacobian->sparsePattern.numberOfNoneZeros; - for(i=0; i < data->simulationInfo->analyticJacobians[index].sizeRows; i++) + for(i=0; i < jacobian->sizeRows; i++) { - data->simulationInfo->analyticJacobians[index].seedVars[i] = 1; + jacobian->seedVars[i] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, systemData->parentJacobian); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - /* infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -data->simulationInfo->analyticJacobians[index].resultVars[l]); */ - systemData->setAElement(i, l, -data->simulationInfo->analyticJacobians[index].resultVars[l], nth, (void*) systemData, threadData); + l = jacobian->sparsePattern.index[ii]; + /* infoStreamPrint(LOG_LS_V, 0, "set on Matrix A (%d, %d)(%d) = %f", i, l, nth, -jacobian->resultVars[l]); */ + systemData->setAElement(i, l, -jacobian->resultVars[l], nth, (void*) systemData, threadData); nth++; ii++; }; @@ -163,7 +165,7 @@ int getAnalyticalJacobianUmfPack(DATA* data, threadData_t *threadData, int sysNu }; /* de-activate seed variable for the corresponding color */ - data->simulationInfo->analyticJacobians[index].seedVars[i] = 0; + jacobian->seedVars[i] = 0; } return 0; @@ -189,7 +191,7 @@ static int wrapper_fvec_umfpack(double* x, double* f, void** data, int sysNumber * author: kbalzereit, wbraun */ int -solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) +solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { void *dataAndThreadData[2] = {data, threadData}; LINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->linearSystemData[sysNumber]); @@ -204,7 +206,6 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) eqSystemNumber, (int) systemData->size, data->localData[0]->timeValue); - rt_ext_tp_tick(&(solverData->timeClock)); if (0 == systemData->method) { @@ -231,7 +232,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) } /* calculate vector b (rhs) */ - memcpy(solverData->work, systemData->x, sizeof(double)*solverData->n_row); + memcpy(solverData->work, aux_x, sizeof(double)*solverData->n_row); wrapper_fvec_umfpack(solverData->work, systemData->b, dataAndThreadData, sysNumber); } tmpJacEvalTime = rt_ext_tp_tock(&(solverData->timeClock)); @@ -242,7 +243,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) { infoStreamPrint(LOG_LS_V, 1, "Old solution x:"); for(i = 0; i < solverData->n_row; ++i) - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); messageClose(LOG_LS_V); infoStreamPrint(LOG_LS_V, 1, "Matrix A n_rows = %d", solverData->n_row); @@ -277,9 +278,9 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) if (0 == status){ if (1 == systemData->method){ - status = umfpack_di_wsolve(UMFPACK_A, solverData->Ap, solverData->Ai, solverData->Ax, systemData->x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); + status = umfpack_di_wsolve(UMFPACK_A, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); } else { - status = umfpack_di_wsolve(UMFPACK_Aat, solverData->Ap, solverData->Ai, solverData->Ax, systemData->x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); + status = umfpack_di_wsolve(UMFPACK_Aat, solverData->Ap, solverData->Ai, solverData->Ax, aux_x, systemData->b, solverData->numeric, solverData->control, solverData->info, solverData->Wi, solverData->W); } } @@ -288,7 +289,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) } else if ((status == UMFPACK_WARNING_singular_matrix) && (casualTearingSet==0)) { - if (!solveSingularSystem(systemData)) + if (!solveSingularSystem(systemData, aux_x)) { success = 1; } @@ -300,10 +301,10 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) if (1 == systemData->method){ /* take the solution */ for(i = 0; i < solverData->n_row; ++i) - systemData->x[i] += solverData->work[i]; + aux_x[i] += solverData->work[i]; /* update inner equations */ - wrapper_fvec_umfpack(systemData->x, solverData->work, dataAndThreadData, sysNumber); + wrapper_fvec_umfpack(aux_x, solverData->work, dataAndThreadData, sysNumber); } else { /* the solution is automatically in x */ } @@ -314,7 +315,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) infoStreamPrint(LOG_LS_V, 0, "System %d numVars %d.", eqSystemNumber, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).numVar); for(i = 0; i < systemData->size; ++i) - infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], systemData->x[i]); + infoStreamPrint(LOG_LS_V, 0, "[%d] %s = %g", i+1, modelInfoGetEquation(&data->modelData->modelDataXml,eqSystemNumber).vars[i], aux_x[i]); messageClose(LOG_LS_V); } @@ -353,7 +354,7 @@ solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber) * * author: kbalzereit, wbraun */ -int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData) +int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData, double* aux_x) { DATA_UMFPACK* solverData = (DATA_UMFPACK*) systemData->solverData[0]; @@ -542,7 +543,7 @@ int solveSingularSystem(LINEAR_SYSTEM_DATA* systemData) /* x = Q^T * z */ for (i = 0; i < solverData->n_col; i++) { - systemData->x[Q[i]] = z[i]; + aux_x[Q[i]] = z[i]; } /* free all used memory */ diff --git a/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.h b/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.h index 92f4e9647f8..7b8d9edc539 100644 --- a/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.h +++ b/SimulationRuntime/c/simulation/solver/linearSolverUmfpack.h @@ -68,7 +68,7 @@ typedef struct DATA_UMFPACK int allocateUmfPackData(int n_row, int n_col, int nz, void **data); int freeUmfPackData(void **data); -int solveUmfPack(DATA *data, threadData_t *threadData, int sysNumber); +int solveUmfPack(DATA *data, threadData_t *threadData, int sysNumbe, double* aux_x); #endif #endif diff --git a/SimulationRuntime/c/simulation/solver/linearSystem.c b/SimulationRuntime/c/simulation/solver/linearSystem.c index d4222bb31b3..66c0734d68e 100644 --- a/SimulationRuntime/c/simulation/solver/linearSystem.c +++ b/SimulationRuntime/c/simulation/solver/linearSystem.c @@ -90,22 +90,22 @@ int initializeLinearSystems(DATA *data, threadData_t *threadData) linsys[i].failed = 0; /* allocate system data */ - linsys[i].x = (double*) malloc(size*sizeof(double)); linsys[i].b = (double*) malloc(size*sizeof(double)); /* check if analytical jacobian is created */ if (1 == linsys[i].method) { + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[linsys[i].jacobianIndex]); if(linsys[i].jacobianIndex != -1) { assertStreamPrint(threadData, 0 != linsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" ); } - if(linsys[i].initialAnalyticalJacobian(data, threadData)) + if(linsys[i].initialAnalyticalJacobian(data, threadData, jacobian)) { linsys[i].jacobianIndex = -1; throwStreamPrint(threadData, "Failed to initialize the jacobian for torn linear system %d.", (int)linsys[i].equationIndex); } - nnz = data->simulationInfo->analyticJacobians[linsys[i].jacobianIndex].sparsePattern.numberOfNoneZeros; + nnz = jacobian->sparsePattern.numberOfNoneZeros; linsys[i].nnz = nnz; } @@ -294,7 +294,6 @@ int freeLinearSystems(DATA *data, threadData_t *threadData) for(i=0; imodelData->nLinearSystems; ++i) { /* free system and solver data */ - free(linsys[i].x); free(linsys[i].b); free(linsys[i].nominal); free(linsys[i].min); @@ -394,7 +393,7 @@ int freeLinearSystems(DATA *data, threadData_t *threadData) * * \author wbraun */ -int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber) +int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x) { TRACE_PUSH int success; @@ -413,7 +412,7 @@ int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber) { #if !defined(OMC_MINIMAL_RUNTIME) case LSS_LIS: - success = solveLis(data, threadData, sysNumber); + success = solveLis(data, threadData, sysNumber, aux_x); break; #else case LSS_LIS_NOT_AVAILABLE: @@ -422,10 +421,10 @@ int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber) #endif #ifdef WITH_UMFPACK case LSS_KLU: - success = solveKlu(data, threadData, sysNumber); + success = solveKlu(data, threadData, sysNumber, aux_x); break; case LSS_UMFPACK: - success = solveUmfPack(data, threadData, sysNumber); + success = solveUmfPack(data, threadData, sysNumber, aux_x); if (!success && linsys->strictTearingFunctionCall != NULL){ debugString(LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used."); success = linsys->strictTearingFunctionCall(data, threadData); @@ -447,20 +446,20 @@ int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber) switch(data->simulationInfo->lsMethod) { case LS_LAPACK: - success = solveLapack(data, threadData, sysNumber); + success = solveLapack(data, threadData, sysNumber, aux_x); break; #if !defined(OMC_MINIMAL_RUNTIME) case LS_LIS: - success = solveLis(data, threadData, sysNumber); + success = solveLis(data, threadData, sysNumber, aux_x); break; #endif #ifdef WITH_UMFPACK case LS_KLU: - success = solveKlu(data, threadData, sysNumber); + success = solveKlu(data, threadData, sysNumber, aux_x); break; case LS_UMFPACK: - success = solveUmfPack(data, threadData, sysNumber); + success = solveUmfPack(data, threadData, sysNumber, aux_x); if (!success && linsys->strictTearingFunctionCall != NULL){ debugString(LOG_DT, "Solving the casual tearing set failed! Now the strict tearing set is used."); success = linsys->strictTearingFunctionCall(data, threadData); @@ -474,11 +473,11 @@ int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber) #endif case LS_TOTALPIVOT: - success = solveTotalPivot(data, threadData, sysNumber); + success = solveTotalPivot(data, threadData, sysNumber, aux_x); break; case LS_DEFAULT: - success = solveLapack(data, threadData, sysNumber); + success = solveLapack(data, threadData, sysNumber, aux_x); /* check if solution process was successful, if not use alternative tearing set if available (dynamic tearing)*/ if (!success && linsys->strictTearingFunctionCall != NULL){ @@ -501,7 +500,7 @@ int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber) logLevel = LOG_STDOUT; } warningStreamPrint(logLevel, 0, "The default linear solver fails, the fallback solver with total pivoting is started at time %f. That might raise performance issues, for more information use -lv LOG_LS.", data->localData[0]->timeValue); - success = solveTotalPivot(data, threadData, sysNumber); + success = solveTotalPivot(data, threadData, sysNumber, aux_x); linsys->failed = 1; }else{ linsys->failed = 0; diff --git a/SimulationRuntime/c/simulation/solver/linearSystem.h b/SimulationRuntime/c/simulation/solver/linearSystem.h index 1e7d0bdcb60..860d6ab3631 100644 --- a/SimulationRuntime/c/simulation/solver/linearSystem.h +++ b/SimulationRuntime/c/simulation/solver/linearSystem.h @@ -51,7 +51,7 @@ typedef void* LS_SOLVER_DATA; int initializeLinearSystems(DATA *data, threadData_t *threadData); int updateStaticDataOfLinearSystems(DATA *data, threadData_t *threadData); int freeLinearSystems(DATA *data, threadData_t *threadData); -int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber); +int solve_linear_system(DATA *data, threadData_t *threadData, int sysNumber, double* aux_x); int check_linear_solutions(DATA *data, int printFailingSystems); void printLinearSystemSolvingStatistics(DATA *data, int sysNumber, int logLevel); diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c b/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c index 53fecfc7694..ce202efe792 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSolverHomotopy.c @@ -807,35 +807,36 @@ int getAnalyticalJacobianHomotopy(DATA_HOMOTOPY* solverData, double* jac) int i,j,k,l,ii; NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[solverData->sysNumber]); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); memset(jac, 0, (solverData->n)*(solverData->n)*sizeof(double)); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { /* activate seed variable for the corresponding color */ - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, NULL); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l; + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; /* Calculate scaled difference quotient */ - jac[k] = data->simulationInfo->analyticJacobians[index].resultVars[l] * solverData->xScaling[j]; + jac[k] = jacobian->resultVars[l] * solverData->xScaling[j]; ii++; }; } /* de-activate seed variable for the corresponding color */ - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[j]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[j] = 0; + if(jacobian->sparsePattern.colorCols[j]-1 == i) + jacobian->seedVars[j] = 0; } } diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c b/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c index 2da098a4446..f7c37eb7471 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSolverHybrd.c @@ -263,35 +263,36 @@ static int getAnalyticalJacobian(struct dataAndSys* dataSys, double* jac) NONLINEAR_SYSTEM_DATA* systemData = &(data->simulationInfo->nonlinearSystemData[dataSys->sysNumber]); DATA_HYBRD* solverData = (DATA_HYBRD*)(systemData->solverData); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); memset(jac, 0, (solverData->n)*(solverData->n)*sizeof(double)); memset(solverData->fjacobian, 0, (solverData->n)*(solverData->n)*sizeof(double)); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { /* activate seed variable for the corresponding color */ - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; - ((systemData->analyticalJacobianColumn))(data, threadData); + ((systemData->analyticalJacobianColumn))(data, threadData, jacobian, NULL); - for(j=0; jsimulationInfo->analyticJacobians[index].sizeCols; j++) + for(j=0; jsizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l; - solverData->fjacobian[k] = jac[k] = data->simulationInfo->analyticJacobians[index].resultVars[l]; + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; + solverData->fjacobian[k] = jac[k] = jacobian->resultVars[l]; ii++; }; } /* de-activate seed variable for the corresponding color */ - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[j]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[j] = 0; + if(jacobian->sparsePattern.colorCols[j]-1 == i) + jacobian->seedVars[j] = 0; } } diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c b/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c index ba079ad4840..94738048373 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSolverNewton.c @@ -82,34 +82,35 @@ int getAnalyticalJacobianNewton(DATA* data, threadData_t *threadData, double* ja NONLINEAR_SYSTEM_DATA* systemData = &(((DATA*)data)->simulationInfo->nonlinearSystemData[currentSys]); DATA_NEWTON* solverData = (DATA_NEWTON*)(systemData->solverData); const int index = systemData->jacobianIndex; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[systemData->jacobianIndex]); memset(jac, 0, (solverData->n)*(solverData->n)*sizeof(double)); - for(i=0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { /* activate seed variable for the corresponding color */ - for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; - systemData->analyticalJacobianColumn(data, threadData); + systemData->analyticalJacobianColumn(data, threadData, jacobian, NULL); - for(j = 0; j < data->simulationInfo->analyticJacobians[index].sizeCols; j++) + for(j = 0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[index].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j]; - while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[j+1]) + ii = jacobian->sparsePattern.leadindex[j]; + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[index].sizeRows + l; - jac[k] = data->simulationInfo->analyticJacobians[index].resultVars[l]; + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; + jac[k] = jacobian->resultVars[l]; ii++; }; } /* de-activate seed variable for the corresponding color */ - if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[j]-1 == i) - data->simulationInfo->analyticJacobians[index].seedVars[j] = 0; + if(jacobian->sparsePattern.colorCols[j]-1 == i) + jacobian->seedVars[j] = 0; } } diff --git a/SimulationRuntime/c/simulation/solver/nonlinearSystem.c b/SimulationRuntime/c/simulation/solver/nonlinearSystem.c index 0462c6f9b30..041ada636e8 100644 --- a/SimulationRuntime/c/simulation/solver/nonlinearSystem.c +++ b/SimulationRuntime/c/simulation/solver/nonlinearSystem.c @@ -380,8 +380,9 @@ int initializeNonlinearSystems(DATA *data, threadData_t *threadData) /* check if analytical jacobian is created */ if(nonlinsys[i].jacobianIndex != -1) { + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[nonlinsys[i].jacobianIndex]); assertStreamPrint(threadData, 0 != nonlinsys[i].analyticalJacobianColumn, "jacobian function pointer is invalid" ); - if(nonlinsys[i].initialAnalyticalJacobian(data, threadData)) + if(nonlinsys[i].initialAnalyticalJacobian(data, threadData, jacobian)) { nonlinsys[i].jacobianIndex = -1; } diff --git a/SimulationRuntime/c/simulation/solver/perform_qss_simulation.c b/SimulationRuntime/c/simulation/solver/perform_qss_simulation.c index ff0658e431a..4f9a201aba0 100644 --- a/SimulationRuntime/c/simulation/solver/perform_qss_simulation.c +++ b/SimulationRuntime/c/simulation/solver/perform_qss_simulation.c @@ -90,12 +90,14 @@ int prefixedName_performQSSSimulation(DATA* data, threadData_t *threadData, SOLV modelica_real *qik, *xik, *derXik, *tq, *tx, *tqp, *nQh, *dQ; modelica_real diffQ = 0.0, dTnextQ = 0.0, nextQ = 0.0; modelica_integer* der = NULL; + const int index = data->callback->INDEX_JAC_A; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[index]); solverInfo->currentTime = simInfo->startTime; warningStreamPrint(LOG_STDOUT, 0, "This QSS method is under development and should not be used yet."); - if (data->callback->initialAnalyticJacobianA(data, threadData)) + if (data->callback->initialAnalyticJacobianA(data, threadData, jacobian)) { infoStreamPrint(LOG_STDOUT, 0, "Jacobian or sparse pattern is not generated or failed to initialize."); return UNKNOWN; diff --git a/SimulationRuntime/c/simulation/solver/stateset.c b/SimulationRuntime/c/simulation/solver/stateset.c index 1517a2352c3..95579ee8d37 100644 --- a/SimulationRuntime/c/simulation/solver/stateset.c +++ b/SimulationRuntime/c/simulation/solver/stateset.c @@ -46,12 +46,17 @@ void initializeStateSetJacobians(DATA *data, threadData_t *threadData) TRACE_PUSH long i = 0; STATE_SET_DATA *set = NULL; + unsigned int jacIndex; + ANALYTIC_JACOBIAN* jacobian; /* go troug all state sets*/ for(i=0; imodelData->nStateSets; i++) { set = &(data->simulationInfo->stateSetData[i]); - if(set->initialAnalyticalJacobian(data, threadData)) + jacIndex = set->jacobianIndex; + jacobian = &(data->simulationInfo->analyticJacobians[jacIndex]); + + if(set->initialAnalyticalJacobian(data, threadData, jacobian)) { throwStreamPrint(threadData, "can not initialze Jacobians for dynamic state selection"); } @@ -139,64 +144,67 @@ static void getAnalyticalJacobianSet(DATA* data, threadData_t *threadData, unsig TRACE_PUSH unsigned int i, j, k, l, ii; unsigned int jacIndex = data->simulationInfo->stateSetData[index].jacobianIndex; - unsigned int nrows = data->simulationInfo->analyticJacobians[jacIndex].sizeRows; - unsigned int ncols = data->simulationInfo->analyticJacobians[jacIndex].sizeCols; + ANALYTIC_JACOBIAN* jacobian = &(data->simulationInfo->analyticJacobians[jacIndex]); + + unsigned int nrows = jacobian->sizeRows; + unsigned int ncols = jacobian->sizeCols; double* jac = data->simulationInfo->stateSetData[index].J; + /* set all elements to zero */ memset(jac, 0, (nrows*ncols*sizeof(double))); - for(i=0; i < data->simulationInfo->analyticJacobians[jacIndex].sparsePattern.maxColors; i++) + for(i=0; i < jacobian->sparsePattern.maxColors; i++) { - for(ii=0; ii < data->simulationInfo->analyticJacobians[jacIndex].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[jacIndex].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[jacIndex].seedVars[ii] = 1; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 1; /* if(ACTIVE_STREAM(LOG_DSS_JAC)) { infoStreamPrint(LOG_DSS_JAC, 1, "Caluculate one col:"); - for(l=0; l < data->simulationInfo->analyticJacobians[jacIndex].sizeCols; l++) - infoStreamPrint(LOG_DSS_JAC, 0, "seed: data->simulationInfo->analyticJacobians[index].seedVars[%d]= %f", l, data->simulationInfo->analyticJacobians[jacIndex].seedVars[l]); + for(l=0; l < jacobian->sizeCols; l++) + infoStreamPrint(LOG_DSS_JAC, 0, "seed: data->simulationInfo->analyticJacobians[index].seedVars[%d]= %f", l, jacobian->seedVars[l]); messageClose(LOG_DSS_JAC); } */ - (data->simulationInfo->stateSetData[index].analyticalJacobianColumn)(data, threadData); + (data->simulationInfo->stateSetData[index].analyticalJacobianColumn)(data, threadData, jacobian, NULL); - for(j=0; j < data->simulationInfo->analyticJacobians[jacIndex].sizeCols; j++) + for(j=0; j < jacobian->sizeCols; j++) { - if(data->simulationInfo->analyticJacobians[jacIndex].seedVars[j] == 1) + if(jacobian->seedVars[j] == 1) { - ii = data->simulationInfo->analyticJacobians[jacIndex].sparsePattern.leadindex[j]; + ii = jacobian->sparsePattern.leadindex[j]; /* infoStreamPrint(LOG_DSS_JAC, 0, "take for %d -> %d\n", j, ii); */ - while(ii < data->simulationInfo->analyticJacobians[jacIndex].sparsePattern.leadindex[j+1]) + while(ii < jacobian->sparsePattern.leadindex[j+1]) { - l = data->simulationInfo->analyticJacobians[jacIndex].sparsePattern.index[ii]; - k = j*data->simulationInfo->analyticJacobians[jacIndex].sizeRows + l; - jac[k] = data->simulationInfo->analyticJacobians[jacIndex].resultVars[l]; - /* infoStreamPrint(LOG_DSS_JAC, 0, "write %d. in jac[%d]-[%d, %d]=%f from col[%d]=%f", ii, k, l, j, jac[k], l, data->simulationInfo->analyticJacobians[jacIndex].resultVars[l]); */ + l = jacobian->sparsePattern.index[ii]; + k = j*jacobian->sizeRows + l; + jac[k] = jacobian->resultVars[l]; + /* infoStreamPrint(LOG_DSS_JAC, 0, "write %d. in jac[%d]-[%d, %d]=%f from col[%d]=%f", ii, k, l, j, jac[k], l, jacobian->resultVars[l]); */ ii++; }; } } - for(ii=0; ii < data->simulationInfo->analyticJacobians[jacIndex].sizeCols; ii++) - if(data->simulationInfo->analyticJacobians[jacIndex].sparsePattern.colorCols[ii]-1 == i) - data->simulationInfo->analyticJacobians[jacIndex].seedVars[ii] = 0; + for(ii=0; ii < jacobian->sizeCols; ii++) + if(jacobian->sparsePattern.colorCols[ii]-1 == i) + jacobian->seedVars[ii] = 0; } if(ACTIVE_STREAM(LOG_DSS_JAC)) { - char *buffer = (char*)malloc(sizeof(char)*data->simulationInfo->analyticJacobians[jacIndex].sizeCols*20); + char *buffer = (char*)malloc(sizeof(char)*jacobian->sizeCols*20); - infoStreamPrint(LOG_DSS_JAC, 1, "jacobian %dx%d [id: %d]", data->simulationInfo->analyticJacobians[jacIndex].sizeRows, data->simulationInfo->analyticJacobians[jacIndex].sizeCols, jacIndex); + infoStreamPrint(LOG_DSS_JAC, 1, "jacobian %dx%d [id: %d]", jacobian->sizeRows, jacobian->sizeCols, jacIndex); - for(i=0; isimulationInfo->analyticJacobians[jacIndex].sizeRows; i++) + for(i=0; isizeRows; i++) { buffer[0] = 0; - for(j=0; j < data->simulationInfo->analyticJacobians[jacIndex].sizeCols; j++) - sprintf(buffer, "%s%.5e ", buffer, jac[i*data->simulationInfo->analyticJacobians[jacIndex].sizeCols+j]); + for(j=0; j < jacobian->sizeCols; j++) + sprintf(buffer, "%s%.5e ", buffer, jac[i*jacobian->sizeCols+j]); infoStreamPrint(LOG_DSS_JAC, 0, "%s", buffer); } messageClose(LOG_DSS_JAC); diff --git a/SimulationRuntime/c/simulation_data.h b/SimulationRuntime/c/simulation_data.h index 6046b769786..da1ea244213 100644 --- a/SimulationRuntime/c/simulation_data.h +++ b/SimulationRuntime/c/simulation_data.h @@ -152,7 +152,6 @@ typedef struct ANALYTIC_JACOBIAN modelica_real* seedVars; modelica_real* tmpVars; modelica_real* resultVars; - modelica_real* jacobian; }ANALYTIC_JACOBIAN; /* EXTERNAL_INPUT @@ -269,8 +268,8 @@ typedef struct NONLINEAR_SYSTEM_DATA * * if analyticalJacobianColumn == NULL no analyticalJacobian is available */ - int (*analyticalJacobianColumn)(void*, threadData_t*); - int (*initialAnalyticalJacobian)(void*, threadData_t*); + int (*analyticalJacobianColumn)(void*, threadData_t*, ANALYTIC_JACOBIAN*, ANALYTIC_JACOBIAN* parentJacobian); + int (*initialAnalyticalJacobian)(void*, threadData_t*, ANALYTIC_JACOBIAN*); modelica_integer jacobianIndex; SPARSE_PATTERN sparsePattern; /* sparse pattern if no jacobian is available */ @@ -321,8 +320,8 @@ typedef struct LINEAR_SYSTEM_DATA void (*setAElement)(int row, int col, double value, int nth, void *data, threadData_t *threadData); void (*setBElement)(int row, double value, void *data, threadData_t *threadData); - int (*analyticalJacobianColumn)(void*, threadData_t*); - int (*initialAnalyticalJacobian)(void*, threadData_t*); + int (*analyticalJacobianColumn)(void*, threadData_t*, ANALYTIC_JACOBIAN*, ANALYTIC_JACOBIAN* parentJacobian); + int (*initialAnalyticalJacobian)(void*, threadData_t*, ANALYTIC_JACOBIAN*); modelica_integer jacobianIndex; void (*residualFunc)(void**, const double*, double*, const int*); @@ -349,6 +348,7 @@ typedef struct LINEAR_SYSTEM_DATA modelica_boolean solved; /* 1: solved in current step - else not */ modelica_boolean failed; /* 1: failed while last try with lapack - else not */ modelica_boolean useSparseSolver; /* 1: use sparse solver, - else any solver */ + ANALYTIC_JACOBIAN* parentJacobian; /* if != NULL then it's the parent jacobian matrix */ /* statistics */ unsigned long numberOfCall; /* number of solving calls of this system */ @@ -405,8 +405,8 @@ typedef struct STATE_SET_DATA * * if analyticalJacobianColumn == NULL no analyticalJacobian is available */ - int (*analyticalJacobianColumn)(void*, threadData_t*); - int (*initialAnalyticalJacobian)(void*, threadData_t*); + int (*analyticalJacobianColumn)(void*, threadData_t*, ANALYTIC_JACOBIAN*, ANALYTIC_JACOBIAN* parentJacobian); + int (*initialAnalyticalJacobian)(void*, threadData_t*, ANALYTIC_JACOBIAN*); modelica_integer jacobianIndex; }STATE_SET_DATA; #else diff --git a/SimulationRuntime/fmi/export/fmi2/fmu2_model_interface.c b/SimulationRuntime/fmi/export/fmi2/fmu2_model_interface.c index 43d09520dcf..0052aac5c5e 100644 --- a/SimulationRuntime/fmi/export/fmi2/fmu2_model_interface.c +++ b/SimulationRuntime/fmi/export/fmi2/fmu2_model_interface.c @@ -487,9 +487,8 @@ fmi2Component fmi2Instantiate(fmi2String instanceName, fmi2Type fmuType, fmi2Str comp->_has_jacobian = 0; comp->fmiDerJac = NULL; if (comp->fmuData->callback->initialPartialFMIDER != NULL){ - if (! comp->fmuData->callback->initialPartialFMIDER(comp->fmuData, comp->threadData)) { + if (! comp->fmuData->callback->initialPartialFMIDER(comp->fmuData, comp->threadData, comp->fmiDerJac)) { comp->_has_jacobian = 1; - comp->fmiDerJac = &(comp->fmuData->simulationInfo->analyticJacobians[comp->fmuData->callback->INDEX_JAC_FMIDER]); } } @@ -1018,7 +1017,7 @@ fmi2Status fmi2GetDirectionalDerivative(fmi2Component c, * More efficient code could only evaluate the equations needed for the * known variables only */ setThreadData(comp); - fmudata->callback->functionJacFMIDER_column(fmudata, td); + fmudata->callback->functionJacFMIDER_column(fmudata, td, comp->fmiDerJac, NULL); resetThreadData(comp); /* Write the results to dvUnknown array */