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

Example 41 with SolverException

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

the class VCellSBMLSolver method solveVCell.

public File solveVCell(String filePrefix, File outDir, String sbmlFileName, SimSpec testSpec) throws IOException, SolverException, SbmlException {
    try {
        cbit.util.xml.VCLogger logger = new LocalLogger();
        // 
        // Instantiate an SBMLImporter to get the speciesUnitsHash - to compute the conversion factor from VC->SB species units.
        // and import SBML  (sbml->bioModel)
        org.vcell.sbml.vcell.SBMLImporter sbmlImporter = new org.vcell.sbml.vcell.SBMLImporter(sbmlFileName, logger, false);
        BioModel bioModel = sbmlImporter.getBioModel();
        if (bRoundTrip) {
            // Round trip the bioModel (bioModel->sbml->bioModel).
            // export bioModel as sbml and save
            String vcml_sbml = cbit.vcell.xml.XmlHelper.exportSBML(bioModel, 2, 1, 0, false, bioModel.getSimulationContext(0), null);
            // re-import bioModel from exported sbml
            XMLSource vcml_sbml_Src = new XMLSource(vcml_sbml);
            BioModel newBioModel = (BioModel) XmlHelper.importSBML(logger, vcml_sbml_Src, false);
            // have rest of code use the round-tripped biomodel
            bioModel = newBioModel;
        }
        // 
        // select only Application, generate math, and create a single Simulation.
        // 
        SimulationContext simContext = bioModel.getSimulationContext(0);
        MathMapping mathMapping = simContext.createNewMathMapping();
        MathDescription mathDesc = mathMapping.getMathDescription();
        simContext.setMathDescription(mathDesc);
        SimulationVersion simVersion = new SimulationVersion(new KeyValue("100"), "unnamed", null, null, null, null, null, null, null, null);
        Simulation sim = new Simulation(simVersion, mathDesc);
        sim.setName("unnamed");
        // if time factor from SBML is not 1 (i.e., it is not in secs but in minutes or hours), convert endTime to min/hr as : endTime*timeFactor
        // double endTime = testSpec.getEndTime()*timeFactor;
        double endTime = testSpec.getEndTime();
        sim.getSolverTaskDescription().setTimeBounds(new TimeBounds(0, endTime));
        TimeStep timeStep = new TimeStep();
        sim.getSolverTaskDescription().setTimeStep(new TimeStep(timeStep.getMinimumTimeStep(), timeStep.getDefaultTimeStep(), endTime / 10000));
        sim.getSolverTaskDescription().setOutputTimeSpec(new UniformOutputTimeSpec((endTime - 0) / testSpec.getNumTimeSteps()));
        sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(testSpec.getAbsTolerance(), testSpec.getRelTolerance()));
        // sim.getSolverTaskDescription().setErrorTolerance(new ErrorTolerance(1e-10, 1e-12));
        // Generate .idaInput string
        File idaInputFile = new File(outDir, filePrefix + SimDataConstants.IDAINPUT_DATA_EXTENSION);
        PrintWriter idaPW = new java.io.PrintWriter(idaInputFile);
        SimulationJob simJob = new SimulationJob(sim, 0, null);
        SimulationTask simTask = new SimulationTask(simJob, 0);
        IDAFileWriter idaFileWriter = new IDAFileWriter(idaPW, simTask);
        idaFileWriter.write();
        idaPW.close();
        // use the idastandalone solver
        File idaOutputFile = new File(outDir, filePrefix + SimDataConstants.IDA_DATA_EXTENSION);
        // String sundialsSolverExecutable = "C:\\Developer\\Eclipse\\workspace\\VCell 4.8\\SundialsSolverStandalone_NoMessaging.exe";
        String executableName = null;
        try {
            executableName = SolverUtilities.getExes(SolverDescription.IDA)[0].getAbsolutePath();
        } catch (IOException e) {
            throw new RuntimeException("failed to get executable for solver " + SolverDescription.IDA.getDisplayLabel() + ": " + e.getMessage(), e);
        }
        Executable executable = new Executable(new String[] { executableName, idaInputFile.getAbsolutePath(), idaOutputFile.getAbsolutePath() });
        executable.start();
        /*			// Generate .cvodeInput string
			File cvodeFile = new File(outDir,filePrefix+SimDataConstants.CVODEINPUT_DATA_EXTENSION);
			PrintWriter cvodePW = new java.io.PrintWriter(cvodeFile);
			SimulationJob simJob = new SimulationJob(sim, 0, null);
		    CVodeFileWriter cvodeFileWriter = new CVodeFileWriter(cvodePW, simJob);
			cvodeFileWriter.write();
			cvodePW.close();

			// use the cvodeStandalone solver
			File cvodeOutputFile = new File(outDir,filePrefix+SimDataConstants.IDA_DATA_EXTENSION);
			String sundialsSolverExecutable = PropertyLoader.getRequiredProperty(PropertyLoader.sundialsSolverExecutableProperty);
			Executable executable = new Executable(new String[]{sundialsSolverExecutable, cvodeFile.getAbsolutePath(), cvodeOutputFile.getAbsolutePath()});
			executable.start();
*/
        // get the result
        ODESolverResultSet odeSolverResultSet = getODESolverResultSet(simJob, idaOutputFile.getPath());
        // remove CVOde input and output files ??
        idaInputFile.delete();
        idaOutputFile.delete();
        // 
        // print header
        // 
        File outputFile = new File(outDir, "results" + filePrefix + ".csv");
        java.io.PrintStream outputStream = new java.io.PrintStream(new java.io.BufferedOutputStream(new java.io.FileOutputStream(outputFile)));
        outputStream.print("time");
        for (int i = 0; i < testSpec.getVarsList().length; i++) {
            outputStream.print("," + testSpec.getVarsList()[i]);
        }
        outputStream.println();
        // 
        // extract data for time and species
        // 
        double[][] data = new double[testSpec.getVarsList().length + 1][];
        int column = odeSolverResultSet.findColumn("t");
        data[0] = odeSolverResultSet.extractColumn(column);
        int origDataLength = data[0].length;
        for (int i = 0; i < testSpec.getVarsList().length; i++) {
            column = odeSolverResultSet.findColumn(testSpec.getVarsList()[i]);
            if (column == -1) {
                Variable var = simJob.getSimulationSymbolTable().getVariable(testSpec.getVarsList()[i]);
                data[i + 1] = new double[data[0].length];
                if (var instanceof cbit.vcell.math.Constant) {
                    double value = ((cbit.vcell.math.Constant) var).getExpression().evaluateConstant();
                    for (int j = 0; j < data[i + 1].length; j++) {
                        data[i + 1][j] = value;
                    }
                } else {
                    throw new RuntimeException("Did not find " + testSpec.getVarsList()[i] + " in simulation");
                }
            } else {
                data[i + 1] = odeSolverResultSet.extractColumn(column);
            }
        }
        // 
        // for each time, print row
        // 
        int index = 0;
        double[] sampleTimes = new double[testSpec.getNumTimeSteps() + 1];
        for (int i = 0; i <= testSpec.getNumTimeSteps(); i++) {
            sampleTimes[i] = endTime * i / testSpec.getNumTimeSteps();
        }
        Model vcModel = bioModel.getModel();
        ReservedSymbol kMole = vcModel.getKMOLE();
        for (int i = 0; i < sampleTimes.length; i++) {
            // 
            while (true) {
                // 
                if (index == odeSolverResultSet.getRowCount() - 1) {
                    if (data[0][index] == sampleTimes[i]) {
                        break;
                    } else {
                        throw new RuntimeException("sampleTime does not match at last time point");
                    }
                }
                // 
                if (data[0][index + 1] > sampleTimes[i]) {
                    break;
                }
                // 
                // sampleTime must be later in our data list.
                // 
                index++;
            }
            // if data[0][index] == sampleTime no need to interpolate
            if (data[0][index] == sampleTimes[i]) {
                // if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
                // if (timeFactor != 1.0) {
                // outputStream.print(data[0][index]/timeFactor);
                // } else {
                outputStream.print(data[0][index]);
                // }
                for (int j = 0; j < testSpec.getVarsList().length; j++) {
                    // SBMLImporter.SBVCConcentrationUnits spConcUnits = speciesUnitsHash.get(testSpec.getVarsList()[j]);
                    // if (spConcUnits != null) {
                    // VCUnitDefinition sbunits = spConcUnits.getSBConcentrationUnits();
                    // VCUnitDefinition vcunits = spConcUnits.getVCConcentrationUnits();
                    // SBMLUnitParameter unitFactor = SBMLUtils.getConcUnitFactor("spConcParam", vcunits, sbunits, kMole);
                    // outputStream.print("," + data[j + 1][index] * unitFactor.getExpression().evaluateConstant()); 		//earlier, hack unitfactor = 0.000001
                    // earlier, hack unitfactor = 0.000001
                    outputStream.print("," + data[j + 1][index]);
                // }
                }
                // System.out.println("No interpolation needed!");
                outputStream.println();
            } else {
                // if data[0][index] < sampleTime, must interpolate
                double fraction = (sampleTimes[i] - data[0][index]) / (data[0][index + 1] - data[0][index]);
                // if timeFactor is not 1.0, time is not in seconds (mins or hrs); if timeFactor is 60, divide sampleTime/60; if it is 3600, divide sampleTime/3600.
                // if (timeFactor != 1.0) {
                // outputStream.print(sampleTimes[i]/timeFactor);
                // } else {
                outputStream.print(sampleTimes[i]);
                // }
                for (int j = 0; j < testSpec.getVarsList().length; j++) {
                    double interpolatedValue = 0.0;
                    double[] speciesVals = null;
                    double[] times = null;
                    // Currently using 2nd order interpolation
                    if (index == 0) {
                        // can only do 1st order interpolation
                        times = new double[] { data[0][index], data[0][index + 1] };
                        speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1] };
                        interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
                    } else if (index >= 1 && index <= origDataLength - 3) {
                        double val_1 = Math.abs(sampleTimes[i] - data[0][index - 1]);
                        double val_2 = Math.abs(sampleTimes[i] - data[0][index + 2]);
                        if (val_1 < val_2) {
                            times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
                            speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
                        } else {
                            times = new double[] { data[0][index], data[0][index + 1], data[0][index + 2] };
                            speciesVals = new double[] { data[j + 1][index], data[j + 1][index + 1], data[j + 1][index + 2] };
                        }
                        interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
                    } else {
                        times = new double[] { data[0][index - 1], data[0][index], data[0][index + 1] };
                        speciesVals = new double[] { data[j + 1][index - 1], data[j + 1][index], data[j + 1][index + 1] };
                        interpolatedValue = MathTestingUtilities.taylorInterpolation(sampleTimes[i], times, speciesVals);
                    }
                    // // Currently using 1st order interpolation
                    // times = new double[] { data[0][index], data[0][index+1] };
                    // speciesVals = new double[] { data[j+1][index], data[j+1][index+1] };
                    // interpolatedValue = taylorInterpolation(sampleTimes[i], times, speciesVals);
                    // interpolatedValue = interpolatedValue * unitFactor.getExpression().evaluateConstant(); 		//earlier, hack unitfactor = 0.000001
                    // System.out.println("Sample time: " + sampleTimes[i] + ", between time[" + index + "]=" + data[0][index]+" and time["+(index+1)+"]="+(data[0][index+1])+", interpolated = "+interpolatedValue);
                    outputStream.print("," + interpolatedValue);
                }
                outputStream.println();
            }
        }
        outputStream.close();
        return outputFile;
    } catch (Exception e) {
        e.printStackTrace(System.out);
        // File outputFile = new File(outDir,"results" + filePrefix + ".csv");
        throw new SolverException(e.getMessage());
    }
}

Example 42 with SolverException

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

the class AdamsMoultonFiveSolver method prep.

/**
 * Integrate over time step using the forward Euler method (1st order explicit)
 * results must be stored in NumVectors-1 = vector(4);
 *  t is the current time
 *  h is the time step
 */
protected void prep(double t, double h) throws SolverException {
    try {
        double[] oldValues = getValueVector(0);
        double[] newValues = getValueVector(1);
        // 
        // update time
        oldValues[getTimeIndex()] = t;
        // newValues has time t, not t + h...it's a
        // scratch array until the end...
        newValues[getTimeIndex()] = t;
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            newValues[I] = oldValues[I];
        }
        for (int i = 0; i < getStateVariableCount(); i++) {
            k[0][getVariableIndex(i)] = h * evaluate(newValues, i);
        }
        // 
        newValues[getTimeIndex()] = t + 0.5 * h;
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            newValues[I] = oldValues[I] + 0.5 * k[0][I];
        }
        for (int i = 0; i < getStateVariableCount(); i++) {
            k[1][getVariableIndex(i)] = h * evaluate(newValues, i);
        }
        // 
        newValues[getTimeIndex()] = t + 0.5 * h;
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            newValues[I] = oldValues[I] + 0.5 * k[1][I];
        }
        for (int i = 0; i < getStateVariableCount(); i++) {
            k[2][getVariableIndex(i)] = h * evaluate(newValues, i);
        }
        // 
        newValues[getTimeIndex()] = t + h;
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            newValues[I] = oldValues[I] + k[2][I];
        }
        for (int i = 0; i < getStateVariableCount(); i++) {
            k[3][getVariableIndex(i)] = h * evaluate(newValues, i);
        }
        // 
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            newValues[I] = oldValues[I] + (k[0][I] + 2.0 * k[1][I] + 2.0 * k[2][I] + k[3][I]) / 6.0;
        }
    } catch (ExpressionException expressionException) {
        throw new SolverException(expressionException.getMessage());
    }
}

Example 43 with SolverException

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

the class DefaultODESolver method integrate.

/**
 * This method was created by a SmartGuide.
 */
protected void integrate() throws SolverException, UserStopException, IOException {
    try {
        SolverTaskDescription taskDescription = simTask.getSimulation().getSolverTaskDescription();
        double timeStep = taskDescription.getTimeStep().getDefaultTimeStep();
        fieldCurrentTime = taskDescription.getTimeBounds().getStartingTime();
        // 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));
        updateResultSet();
        // 
        int iteration = 0;
        while (fieldCurrentTime < taskDescription.getTimeBounds().getEndingTime()) {
            checkForUserStop();
            step(fieldCurrentTime, timeStep);
            // 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 += timeStep;
            iteration++;
            fieldCurrentTime = taskDescription.getTimeBounds().getStartingTime() + iteration * timeStep;
            // Print results if it coincides with a save interval...
            if (taskDescription.getOutputTimeSpec().isDefault()) {
                int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
                if ((iteration % keepEvery) == 0) {
                    updateResultSet();
                }
            }
        }
        // store last time point
        if (taskDescription.getOutputTimeSpec().isDefault()) {
            int keepEvery = ((DefaultOutputTimeSpec) taskDescription.getOutputTimeSpec()).getKeepEvery();
            if ((iteration % keepEvery) != 0)
                updateResultSet();
        }
    } catch (ExpressionException | MathException e) {
        throw new SolverException("Solver failed: " + e.getMessage(), e);
    }
}

Example 44 with SolverException

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

the class ForwardEulerSolver method step.

/**
 * Integrate over time step using the forward Euler method (1st order explicit)
 * results must be stored in NumVectors-1 = vector(1);
 */
protected void step(double t, double h) throws SolverException {
    try {
        // get value Vectors
        double[] oldValues = getValueVector(0);
        double[] newValues = getValueVector(1);
        // update time
        oldValues[getTimeIndex()] = t;
        newValues[getTimeIndex()] = t + h;
        // integrate
        for (int i = 0; i < getStateVariableCount(); i++) {
            int I = getVariableIndex(i);
            newValues[I] = oldValues[I] + h * evaluate(oldValues, i);
        }
    } catch (ExpressionException expressionException) {
        throw new SolverException(expressionException.getMessage());
    }
}

Example 45 with SolverException

use of cbit.vcell.solver.SolverException 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());
    }
}