Examples with SolverTaskDescription cbit.vcell.solver.SolverTaskDescription used on opensource projects

Example 36 with SolverTaskDescription

use of cbit.vcell.solver.SolverTaskDescription in project vcell by virtualcell.

the class OdeFileWriter method write.

/**
 * Insert the method's description here.
 * Creation date: (3/8/00 10:31:52 PM)
 */
public void write(String[] parameterNames) throws Exception {
    createStateVariables();
    createSymbolTable();
    Simulation simulation = simTask.getSimulation();
    if (simulation.getSolverTaskDescription().getUseSymbolicJacobian()) {
        throw new RuntimeException("symbolic jacobian option not yet supported in interpreted Stiff solver");
    }
    writeJMSParamters();
    SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
    TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
    ErrorTolerance errorTolerance = solverTaskDescription.getErrorTolerance();
    printWriter.println("SOLVER " + getSolverName());
    printWriter.println("STARTING_TIME " + timeBounds.getStartingTime());
    printWriter.println("ENDING_TIME " + timeBounds.getEndingTime());
    printWriter.println("RELATIVE_TOLERANCE " + errorTolerance.getRelativeErrorTolerance());
    printWriter.println("ABSOLUTE_TOLERANCE " + errorTolerance.getAbsoluteErrorTolerance());
    printWriter.println("MAX_TIME_STEP " + simulation.getSolverTaskDescription().getTimeStep().getMaximumTimeStep());
    OutputTimeSpec ots = simulation.getSolverTaskDescription().getOutputTimeSpec();
    if (ots.isDefault()) {
        printWriter.println("KEEP_EVERY " + ((DefaultOutputTimeSpec) ots).getKeepEvery());
    } else if (ots.isUniform()) {
        printWriter.println("OUTPUT_TIME_STEP " + ((UniformOutputTimeSpec) ots).getOutputTimeStep());
    } else if (ots.isExplicit()) {
        printWriter.println("OUTPUT_TIMES " + ((ExplicitOutputTimeSpec) ots).getNumTimePoints());
        printWriter.println(((ExplicitOutputTimeSpec) ots).toSpaceSeperatedMultiLinesOfString());
    }
    if (parameterNames != null && parameterNames.length != 0) {
        printWriter.println("NUM_PARAMETERS " + parameterNames.length);
        for (int i = 0; i < parameterNames.length; i++) {
            printWriter.println(parameterNames[i]);
        }
    }
    HashMap<Discontinuity, String> discontinuityNameMap = new HashMap<Discontinuity, String>();
    String eventString = null;
    if (simulation.getMathDescription().hasEvents()) {
        eventString = writeEvents(discontinuityNameMap);
    }
    String equationString = writeEquations(discontinuityNameMap);
    if (discontinuityNameMap.size() > 0) {
        printWriter.println("DISCONTINUITIES " + discontinuityNameMap.size());
        for (Discontinuity od : discontinuityNameMap.keySet()) {
            printWriter.println(discontinuityNameMap.get(od) + " " + od.getDiscontinuityExp().flatten().infix() + "; " + od.getRootFindingExp().flatten().infix() + ";");
        }
    }
    if (eventString != null) {
        printWriter.print(eventString);
    }
    printWriter.println("NUM_EQUATIONS " + getStateVariableCount());
    printWriter.println(equationString);
}

Example 37 with SolverTaskDescription

use of cbit.vcell.solver.SolverTaskDescription in project vcell by virtualcell.

the class RungeKuttaFehlbergSolver method integrate.

protected void integrate() throws cbit.vcell.solver.SolverException, UserStopException, IOException {
    try {
        // Get machine epsilon
        final double DBL_EPSILON = 1.0E-12;
        final double epsilon = DBL_EPSILON;
        final double twentySixEpsilon = 26 * epsilon;
        // 
        SolverTaskDescription taskDescription = simTask.getSimulation().getSolverTaskDescription();
        double startingTime = taskDescription.getTimeBounds().getStartingTime();
        double endingTime = taskDescription.getTimeBounds().getEndingTime();
        double relativeErrorTolerance = taskDescription.getErrorTolerance().getRelativeErrorTolerance();
        double absoluteErrorTolerance = taskDescription.getErrorTolerance().getAbsoluteErrorTolerance();
        fieldCurrentTime = taskDescription.getTimeBounds().getStartingTime();
        // 
        double[] oldValues = getValueVector(0);
        double[] newValues = getValueVector(1);
        // The minimumRelativeError is the minimum acceptable value of
        // relativeError.  Attempts to obtain higher accuracy with this
        // subroutine are usually very expensive and often unsuccessful.
        final double minimumRelativeError = 1e-12;
        // Check input parameters...
        if (relativeErrorTolerance < 0.0 || absoluteErrorTolerance < 0.0 || startingTime >= endingTime) {
            throw new SolverException("Invalid parameters");
        }
        // difficulties arising from impossible accuracy requests
        if (relativeErrorTolerance < 2 * epsilon + minimumRelativeError) {
            // Or just set relativeError = 2 * epsilon + minimumRelativeError???
            throw new SolverException("Relative error too small");
        }
        // before computation begins, settle fast equilibrium
        if (getFastAlgebraicSystem() != null) {
            fieldValueVectors.copyValues(0, 1);
            getFastAlgebraicSystem().initVars(getValueVector(0), getValueVector(1));
            getFastAlgebraicSystem().solveSystem(getValueVector(0), getValueVector(1));
            fieldValueVectors.copyValues(1, 0);
        }
        // check for failure
        check(getValueVector(0));
        // 
        double timeRemaining = endingTime - fieldCurrentTime;
        double tolerance = 0.0;
        double maximumTolerance = 0.0;
        final double maximumTimeStep = taskDescription.getTimeStep().getMaximumTimeStep();
        double h = Math.min(Math.abs(timeRemaining), taskDescription.getTimeStep().getMaximumTimeStep());
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            tolerance = relativeErrorTolerance * Math.abs(oldValues[I]) + absoluteErrorTolerance;
            if (tolerance > 0.0) {
                double yp = Math.abs(evaluate(oldValues, i));
                if (yp * Math.pow(h, 5.0) > tolerance) {
                    h = Math.pow(tolerance / yp, 0.2);
                }
            }
            maximumTolerance = Math.max(maximumTolerance, tolerance);
        }
        if (maximumTolerance <= 0.0)
            h = 0.0;
        h = Math.max(h, twentySixEpsilon * Math.max(timeRemaining, Math.abs(fieldCurrentTime)));
        // To avoid premature underflow in the error
        // tolerance function,  scale the tolerances...
        final double scale = 2.0 / relativeErrorTolerance;
        final double scaledAbsoluteError = scale * absoluteErrorTolerance;
        boolean previousStepFailed = false;
        // Check for failure...
        check(getValueVector(0));
        updateResultSet();
        // Integrate...
        int iteration = 0;
        while (fieldCurrentTime < endingTime) {
            checkForUserStop();
            // Set smallest allowable stepsize...
            final double minimumTimeStep = Math.max(twentySixEpsilon * Math.abs(fieldCurrentTime), taskDescription.getTimeStep().getMinimumTimeStep());
            h = Math.min(h, taskDescription.getTimeStep().getMaximumTimeStep());
            // Adjust stepsize if necessary to hit the output point.
            // Look ahead two steps to avoid drastic changes in the stepsize and
            // thus lessen the impact of output points on the code.
            timeRemaining = endingTime - fieldCurrentTime;
            if (timeRemaining < 2 * h) {
                if (timeRemaining <= h) {
                    h = timeRemaining;
                } else {
                    h = 0.5 * timeRemaining;
                }
            }
            // If too close to output point, extrapolate and return
            if (timeRemaining <= twentySixEpsilon * Math.abs(fieldCurrentTime)) {
                for (int i = 0; i < getStateVariableCount(); i++) {
                    int I = getVariableIndex(i);
                    newValues[I] = oldValues[I] + timeRemaining * evaluate(oldValues, i);
                }
                // update (old = new)
                fieldValueVectors.copyValuesDown();
                // compute fast system
                if (getFastAlgebraicSystem() != null) {
                    fieldValueVectors.copyValues(0, 1);
                    getFastAlgebraicSystem().initVars(getValueVector(0), getValueVector(1));
                    getFastAlgebraicSystem().solveSystem(getValueVector(0), getValueVector(1));
                    fieldValueVectors.copyValues(1, 0);
                }
                // check for failure
                check(getValueVector(0));
                fieldCurrentTime = endingTime;
                return;
            } else {
                // Advance an approximate solution over one step of length h
                // RungeKuttaStep(Model,y,t,h,ss);
                step(fieldCurrentTime, h);
                // compute fast system
                if (getFastAlgebraicSystem() != null) {
                    fieldValueVectors.copyValues(1, 2);
                    getFastAlgebraicSystem().initVars(getValueVector(1), getValueVector(2));
                    getFastAlgebraicSystem().solveSystem(getValueVector(1), getValueVector(2));
                    fieldValueVectors.copyValues(2, 1);
                }
                // check for failure
                check(getValueVector(1));
                // Compute and test allowable tolerances versus local error estimates
                // and remove scaling of tolerances. Note that relative error is
                // measured with respect to the average of the magnitudes of the
                // solution at the beginning and end of the step.
                double estimatedErrorOverErrorTolerance = 0.0;
                for (int i = 0; i < getStateVariableCount(); i++) {
                    int I = getVariableIndex(i);
                    double errorTolerance = Math.abs(oldValues[I]) + Math.abs(newValues[I]) + scaledAbsoluteError;
                    // Inappropriate error tolerance
                    if (errorTolerance <= 0.0) {
                        throw new SolverException("Error tolerance too small");
                    }
                    double estimatedError = Math.abs(calculateErrorTerm(i));
                    estimatedErrorOverErrorTolerance = Math.max(estimatedErrorOverErrorTolerance, estimatedError / errorTolerance);
                }
                double estimatedTolerance = h * estimatedErrorOverErrorTolerance * scale;
                if (estimatedTolerance > 1.0) {
                    // Unsuccessful step.  Reduce the stepsize and try again.
                    // The decrease is limited to a factor of 1/10...
                    previousStepFailed = true;
                    double s = 0.1;
                    if (estimatedTolerance < 59049.0) {
                        // s = 0.1  @  estimatedTolerance = 59049
                        s = 0.9 / Math.pow(estimatedTolerance, 0.2);
                    }
                    h *= s;
                    if (h < minimumTimeStep) {
                        throw new SolverException("Requested error unattainable at smallest allowable stepsize");
                    }
                } else {
                    // Successful step.  Store solution at t+h and evaluate derivatives there.
                    fieldValueVectors.copyValuesDown();
                    fieldCurrentTime += h;
                    iteration++;
                    if (taskDescription.getOutputTimeSpec().isDefault()) {
                        int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
                        if ((iteration % keepEvery) == 0)
                            updateResultSet();
                    }
                    // 
                    // Choose next stepsize.  The increase is limited to a factor of 5. If
                    // step failure has just occurred, next stepsize is not allowed to increase.
                    // 
                    double s = 5;
                    if (estimatedTolerance > 0.0001889568) {
                        // s = 5  @  estimatedTolerance = 0.0001889568
                        s = 0.9 / Math.pow(estimatedTolerance, 0.2);
                    }
                    if (previousStepFailed) {
                        s = Math.min(1.0, s);
                    }
                    h = Math.min(Math.max(minimumTimeStep, s * h), maximumTimeStep);
                    previousStepFailed = false;
                }
            }
        }
        // store last time point
        if (taskDescription.getOutputTimeSpec().isDefault()) {
            int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
            if ((iteration % keepEvery) != 0)
                updateResultSet();
        }
    } catch (ExpressionException expressionException) {
        expressionException.printStackTrace(System.out);
        throw new SolverException(expressionException.getMessage());
    } catch (MathException mathException) {
        mathException.printStackTrace(System.out);
        throw new SolverException(mathException.getMessage());
    }
}

Example 38 with SolverTaskDescription

use of cbit.vcell.solver.SolverTaskDescription in project vcell by virtualcell.

the class HtcSimulationWorker method choresFor.

/**
 * determine post processing chores to been done after the simulation completes
 * @param simTask
 * @return PostProcessingChores
 */
private PostProcessingChores choresFor(SimulationTask simTask) {
    String userDir = "/" + simTask.getUserName();
    String primaryInternal = PropertyLoader.getRequiredProperty(PropertyLoader.primarySimDataDirInternalProperty);
    String primaryExternal = PropertyLoader.getRequiredProperty(PropertyLoader.primarySimDataDirExternalProperty);
    PostProcessingChores chores = null;
    final SolverTaskDescription slvTaskDesc = simTask.getSimulation().getSolverTaskDescription();
    if (!slvTaskDesc.isParallel()) {
        chores = new PostProcessingChores(primaryInternal + userDir, primaryExternal + userDir);
    } else {
        String runDirExternal = PropertyLoader.getRequiredProperty(PropertyLoader.PARALLEL_DATA_DIR_EXTERNAL);
        chores = new PostProcessingChores(runDirExternal + userDir, primaryExternal + userDir);
    }
    chores.setVtkUser(slvTaskDesc.isVtkUser());
    if (lg.isDebugEnabled()) {
        lg.debug("Simulation " + simTask.getSimulation().getDescription() + " task " + simTask.getTaskID() + " with " + slvTaskDesc.getNumProcessors() + " processors using " + chores);
    }
    return chores;
}

Example 39 with SolverTaskDescription

use of cbit.vcell.solver.SolverTaskDescription in project vcell by virtualcell.

the class SimulationTable method getSQLValueList.

/**
 * This method was created in VisualAge.
 * @return java.lang.String
 * @param key KeyValue
 * @param modelName java.lang.String
 */
public String getSQLValueList(Simulation simulation, KeyValue mathKey, Version version, DatabaseSyntax dbSyntax) {
    SolverTaskDescription solverTD = simulation.getSolverTaskDescription();
    String mathOverrides = simulation.getMathOverrides().getVCML();
    StringBuffer buffer = new StringBuffer();
    switch(dbSyntax) {
        case ORACLE:
            {
                buffer.append("(");
                buffer.append(getVersionGroupSQLValue(version) + ",");
                buffer.append(mathKey + ",");
                // MathOverridesOrig keep for compatibility with old VCell
                buffer.append("EMPTY_CLOB()" + ",");
                break;
            }
        case POSTGRES:
            {
                buffer.append("(");
                buffer.append(getVersionGroupSQLValue(version) + ",");
                buffer.append(mathKey + ",");
                // MathOverridesOrig keep for compatibility with old VCell
                buffer.append("null" + ",");
                break;
            }
        default:
            {
                throw new RuntimeException("unexpected DatabaseSyntax " + dbSyntax);
            }
    }
    if (DbDriver.varchar2_CLOB_is_Varchar2_OK(mathOverrides)) {
        buffer.append("null" + "," + DbDriver.INSERT_VARCHAR2_HERE + ",");
    } else {
        buffer.append(DbDriver.INSERT_CLOB_HERE + "," + "null" + ",");
    }
    buffer.append((solverTD != null ? "'" + TokenMangler.getSQLEscapedString(solverTD.getVCML()) + "'" : "null") + ",");
    if (simulation.getMathDescription() != null && simulation.getMathDescription().getGeometry() != null && simulation.getMathDescription().getGeometry().getDimension() > 0) {
        MeshSpecification meshSpec = simulation.getMeshSpecification();
        buffer.append(meshSpec.getSamplingSize().getX() + "," + meshSpec.getSamplingSize().getY() + "," + meshSpec.getSamplingSize().getZ());
    } else {
        buffer.append("null,null,null");
    }
    if (simulation.getDataProcessingInstructions() != null) {
        buffer.append(",'" + TokenMangler.getSQLEscapedString(simulation.getDataProcessingInstructions().getDbXml()) + "'");
    } else {
        buffer.append(",null");
    }
    buffer.append(")");
    return buffer.toString();
}