Examples with TimeInterval org.vcell.chombo.TimeInterval used on opensource projects

Example 1 with TimeInterval

use of org.vcell.chombo.TimeInterval in project vcell by virtualcell.

the class FiniteVolumeFileWriter method writeSimulationParamters.

/**
 * @param timeFunction
 * @param startTime
 * @param endTime
 * @param rootFinder
 * @param uniqueRootTimes
 * @param bPrintIterations
 * @throws ExpressionException
 *
 * for testing within scrapbook, see below:
 *
 * try {
 *	edu.northwestern.at.utils.math.rootfinders.MonadicFunctionRootFinder rootFinder =
 *//			new edu.northwestern.at.utils.math.rootfinders.Brent();
 *			new edu.northwestern.at.utils.math.rootfinders.Bisection();
 *//			new edu.northwestern.at.utils.math.rootfinders.NewtonRaphson();
 *//			new edu.northwestern.at.utils.math.rootfinders.Secant();
 *	cbit.vcell.parser.SimpleSymbolTable simpleSymbolTable = new cbit.vcell.parser.SimpleSymbolTable(new String[] { "t" });
 *
 *	cbit.vcell.parser.Expression exp = new cbit.vcell.parser.Expression("t-0.56");
 *
 *	exp.bindExpression(simpleSymbolTable);
 *	java.util.TreeSet<Double> rootTimes = new java.util.TreeSet<Double>();
 *	double startTime = 0.0;
 *	double endTime = 100.0;
 *	System.out.print("exp = '"+ exp.infix() + "'");
 *	long currentTimeMS = System.currentTimeMillis();
 *	cbit.vcell.solvers.FiniteVolumeFileWriter.findAllRoots(exp,startTime,endTime,rootFinder,rootTimes,false);
 *	long finalTimeMS = System.currentTimeMillis();
 *	for (double root : rootTimes){
 *		System.out.println("root = "+root);
 *	}
 *	System.out.println("elapsedTime of computation = "+(finalTimeMS-currentTimeMS)+" ms, found " + rootTimes.size() + " roots (not unique)");
 *
 *}catch (Exception e){
 *	e.printStackTrace(System.out);
 *}
 */
/*public static void findAllRoots(Expression timeFunction, double startTime, double endTime, MonadicFunctionRootFinder rootFinder, TreeSet<Double> uniqueRootTimes, boolean bPrintIterations) throws ExpressionException{
	TreeSet<Double> allRootTimes = new TreeSet<Double>();
	final Expression function_exp = new Expression(timeFunction);
	MonadicFunction valueFunction = new MonadicFunction() {
		double[] values = new double[1];
		public double f(double t) {
			values[0] = t;
			try {
				return function_exp.evaluateVector(values);
			} catch (ExpressionException e) {
				e.printStackTrace();
				throw new RuntimeException("expression exception "+e.getMessage());
			}
		}
	};
	
	final Expression derivative_exp = new Expression(timeFunction.differentiate(ReservedVariable.TIME.getName()));
	MonadicFunction derivativeFunction = new MonadicFunction() {
		double[] values = new double[1];
		public double f(double t) {
			values[0] = t;
			try {
				return derivative_exp.evaluateVector(values);
			} catch (ExpressionException e) {
				e.printStackTrace();
				throw new RuntimeException("expression exception "+e.getMessage());
			}
		}
	};
	
	RootFinderConvergenceTest convergenceTest = new StandardRootFinderConvergenceTest();
	RootFinderIterationInformation iterationInformation = null;
	if (bPrintIterations){
		iterationInformation = new RootFinderIterationInformation() {				
			public void iterationInformation(double x, double fx, double dfx, int currentIteration) {
				System.out.println(currentIteration+") x="+x+", fx="+fx+", dfx="+dfx);
			}
		};
	}
	int NUM_BRACKETS = 1000;
	double simulationTime = endTime - startTime;
	double tolerance = simulationTime/1e10;
	int maxIter = 1000;
	
	for (int i=0;i<NUM_BRACKETS-1;i++){
		double bracketMin = startTime + simulationTime*i/NUM_BRACKETS;
		double bracketMax = startTime + simulationTime*(i+1)/NUM_BRACKETS;
	
		double root = rootFinder.findRoot(bracketMin, bracketMax, tolerance, maxIter, valueFunction, derivativeFunction, convergenceTest, iterationInformation);
		if (root>startTime && root<endTime && valueFunction.f(root)<=tolerance){
			allRootTimes.add(root);
		}
	}
	double uniqueTolerance = tolerance * 100;
	double lastUniqueRoot = Double.NEGATIVE_INFINITY;
	for (double root : allRootTimes){
		if (root-lastUniqueRoot > uniqueTolerance){
			uniqueRootTimes.add(root);
		}
		lastUniqueRoot = root;
	}

}  ---------------------------------JIM's CODE COMMENTTED FOR FUTURE DEVELOPMENT*/
/**
 *# Simulation Parameters
 *SIMULATION_PARAM_BEGIN
 *SOLVER SUNDIALS_PDE_SOLVER 1.0E-7 1.0E-9
 *DISCONTINUITY_TIMES 2 1.0E-4 3.0000000000000003E-4
 *BASE_FILE_NAME c:/Vcell/users/fgao/SimID_31746636_0_
 *ENDING_TIME 4.0E-4
 *KEEP_EVERY ONE_STEP 3
 *KEEP_AT_MOST 1000
 *SIMULATION_PARAM_END
 * @throws MathException
 * @throws ExpressionException
 */
private void writeSimulationParamters() throws ExpressionException, MathException {
    Simulation simulation = simTask.getSimulation();
    SolverTaskDescription solverTaskDesc = simulation.getSolverTaskDescription();
    printWriter.println("# Simulation Parameters");
    printWriter.println(FVInputFileKeyword.SIMULATION_PARAM_BEGIN);
    if (solverTaskDesc.getSolverDescription().equals(SolverDescription.SundialsPDE)) {
        printWriter.print(FVInputFileKeyword.SOLVER + " " + FVInputFileKeyword.SUNDIALS_PDE_SOLVER + " " + solverTaskDesc.getErrorTolerance().getRelativeErrorTolerance() + " " + solverTaskDesc.getErrorTolerance().getAbsoluteErrorTolerance() + " " + solverTaskDesc.getTimeStep().getMaximumTimeStep());
        if (simulation.getMathDescription().hasVelocity()) {
            printWriter.print(" " + solverTaskDesc.getSundialsPdeSolverOptions().getMaxOrderAdvection());
        }
        printWriter.println();
        Vector<Discontinuity> discontinuities = new Vector<Discontinuity>();
        TreeSet<Double> discontinuityTimes = new TreeSet<Double>();
        MathDescription mathDesc = simulation.getMathDescription();
        Enumeration<SubDomain> enum1 = mathDesc.getSubDomains();
        while (enum1.hasMoreElements()) {
            SubDomain sd = enum1.nextElement();
            Enumeration<Equation> enum_equ = sd.getEquations();
            while (enum_equ.hasMoreElements()) {
                Equation equation = enum_equ.nextElement();
                equation.getDiscontinuities(simTask.getSimulationJob().getSimulationSymbolTable(), discontinuities);
            }
        }
        getDiscontinuityTimes(discontinuities, discontinuityTimes);
        if (discontinuityTimes.size() > 0) {
            printWriter.print(FVInputFileKeyword.DISCONTINUITY_TIMES + " " + discontinuityTimes.size());
            for (double d : discontinuityTimes) {
                printWriter.print(" " + d);
            }
            printWriter.println();
        }
    } else if (solverTaskDesc.getSolverDescription().equals(SolverDescription.Chombo)) {
        printWriter.println(FVInputFileKeyword.SOLVER + " " + FVInputFileKeyword.CHOMBO_SEMIIMPLICIT_SOLVER);
    } else if (solverTaskDesc.getSolverDescription().equals(SolverDescription.VCellPetsc)) {
        printWriter.println(FVInputFileKeyword.SOLVER + " " + FVInputFileKeyword.VCELL_PETSC_SOLVER);
    } else {
        printWriter.println(FVInputFileKeyword.SOLVER + " " + FVInputFileKeyword.FV_SOLVER + " " + solverTaskDesc.getErrorTolerance().getRelativeErrorTolerance());
    }
    printWriter.println(FVInputFileKeyword.BASE_FILE_NAME + " " + new File(workingDirectory, simTask.getSimulationJob().getSimulationJobID()).getAbsolutePath());
    if (solverTaskDesc.isParallel() && destinationDirectory != null && !destinationDirectory.equals(workingDirectory)) {
        printWriter.println(FVInputFileKeyword.PRIMARY_DATA_DIR + " " + destinationDirectory.getAbsolutePath());
    }
    printWriter.println(FVInputFileKeyword.ENDING_TIME + " " + solverTaskDesc.getTimeBounds().getEndingTime());
    OutputTimeSpec outputTimeSpec = solverTaskDesc.getOutputTimeSpec();
    if (solverTaskDesc.getSolverDescription().isChomboSolver()) {
        List<TimeInterval> timeIntervalList = solverTaskDesc.getChomboSolverSpec().getTimeIntervalList();
        printWriter.println(FVInputFileKeyword.TIME_INTERVALS + " " + timeIntervalList.size());
        for (TimeInterval ti : timeIntervalList) {
            printWriter.println(ti.getEndingTime() + " " + ti.getTimeStep() + " " + ti.getKeepEvery());
        }
    } else if (solverTaskDesc.getSolverDescription().equals(SolverDescription.SundialsPDE)) {
        if (outputTimeSpec.isDefault()) {
            DefaultOutputTimeSpec defaultOutputTimeSpec = (DefaultOutputTimeSpec) outputTimeSpec;
            printWriter.println(FVInputFileKeyword.KEEP_EVERY + " " + FVInputFileKeyword.ONE_STEP + " " + defaultOutputTimeSpec.getKeepEvery());
            printWriter.println(FVInputFileKeyword.KEEP_AT_MOST + " " + defaultOutputTimeSpec.getKeepAtMost());
        } else {
            printWriter.println(FVInputFileKeyword.TIME_STEP + " " + ((UniformOutputTimeSpec) outputTimeSpec).getOutputTimeStep());
            printWriter.println(FVInputFileKeyword.KEEP_EVERY + " 1");
        }
    } else {
        double defaultTimeStep = solverTaskDesc.getTimeStep().getDefaultTimeStep();
        printWriter.println(FVInputFileKeyword.TIME_STEP + " " + defaultTimeStep);
        int keepEvery = 1;
        if (outputTimeSpec.isDefault()) {
            keepEvery = ((DefaultOutputTimeSpec) outputTimeSpec).getKeepEvery();
        } else if (outputTimeSpec.isUniform()) {
            UniformOutputTimeSpec uots = (UniformOutputTimeSpec) outputTimeSpec;
            double ots = uots.getOutputTimeStep();
            keepEvery = (int) Math.round(ots / defaultTimeStep);
        } else {
            throw new RuntimeException("unexpected OutputTime specification type :" + outputTimeSpec.getClass().getName());
        }
        if (keepEvery <= 0) {
            throw new RuntimeException("Output KeepEvery must be a positive integer. Try to change the output option.");
        }
        printWriter.println(FVInputFileKeyword.KEEP_EVERY + " " + keepEvery);
    }
    ErrorTolerance stopAtSpatiallyUniformErrorTolerance = solverTaskDesc.getStopAtSpatiallyUniformErrorTolerance();
    if (stopAtSpatiallyUniformErrorTolerance != null) {
        printWriter.println(FVInputFileKeyword.CHECK_SPATIALLY_UNIFORM + " " + stopAtSpatiallyUniformErrorTolerance.getAbsoluteErrorTolerance() + " " + stopAtSpatiallyUniformErrorTolerance.getRelativeErrorTolerance());
    }
    printWriter.println(FVInputFileKeyword.SIMULATION_PARAM_END);
    printWriter.println();
}

Example 2 with TimeInterval

use of org.vcell.chombo.TimeInterval in project vcell by virtualcell.

the class XmlReader method getChomboSolverSpec.

private ChomboSolverSpec getChomboSolverSpec(SolverTaskDescription solverTaskDesc, Element element, int dimension) throws XmlParseException {
    int maxBoxSize = parseIntWithDefault(element, XMLTags.MaxBoxSizeTag, ChomboSolverSpec.getDefaultMaxBoxSize(dimension));
    double fillRatio = parseDoubleWithDefault(element, XMLTags.FillRatioTag, ChomboSolverSpec.getDefaultFillRatio());
    boolean bSaveVCellOutput = parseBooleanWithDefault(element, XMLTags.SaveVCellOutput, true);
    boolean bSaveChomboOutput = parseBooleanWithDefault(element, XMLTags.SaveChomboOutput, false);
    Element childElement = element.getChild(XMLTags.RefineRatios, vcNamespace);
    List<Integer> refineRatioList = null;
    if (childElement != null) {
        String text = childElement.getText();
        if (text != null && !text.isEmpty()) {
            StringTokenizer st = new StringTokenizer(text, ",");
            if (st.hasMoreTokens()) {
                refineRatioList = new ArrayList<Integer>();
                while (st.hasMoreElements()) {
                    String token = st.nextToken();
                    if (token != null) {
                        int n = Integer.parseInt(token);
                        refineRatioList.add(n);
                    }
                }
            }
        }
    }
    Integer viewLevel = null;
    try {
        ChomboSolverSpec css = new ChomboSolverSpec(maxBoxSize, fillRatio, viewLevel, bSaveVCellOutput, bSaveChomboOutput, refineRatioList);
        double smallVolfracThreshold = parseDoubleWithDefault(element, XMLTags.SmallVolfracThreshold, 0);
        int blockFactor = parseIntWithDefault(element, XMLTags.BlockFactorTag, ChomboSolverSpec.DEFAULT_BLOCK_FACTOR);
        boolean bActivateFeatureUnderDevelopment = parseBooleanWithDefault(element, XMLTags.ActivateFeatureUnderDevelopment, false);
        css.setSmallVolfracThreshold(smallVolfracThreshold);
        css.setActivateFeatureUnderDevelopment(bActivateFeatureUnderDevelopment);
        css.setBlockFactor(blockFactor);
        int tagsGrow = parseIntWithDefault(element, XMLTags.TagsGrowTag, ChomboSolverSpec.defaultTagsGrow);
        css.setTagsGrow(tagsGrow);
        Element timeBoundsElement = element.getChild(XMLTags.TimeBoundTag, vcNamespace);
        List<Element> timeIntervalElementList = null;
        boolean noTimeBounds = false;
        if (timeBoundsElement == null) {
            noTimeBounds = true;
        } else {
            timeIntervalElementList = timeBoundsElement.getChildren(XMLTags.TimeIntervalTag, vcNamespace);
            if (timeIntervalElementList.size() == 0) {
                noTimeBounds = true;
            }
        }
        if (noTimeBounds) {
            // old format
            double startTime = 0;
            double endTime = solverTaskDesc.getTimeBounds().getEndingTime();
            double timeStep = solverTaskDesc.getTimeStep().getDefaultTimeStep();
            double outputTimeStep = ((UniformOutputTimeSpec) solverTaskDesc.getOutputTimeSpec()).getOutputTimeStep();
            try {
                TimeInterval ti = new TimeInterval(startTime, endTime, timeStep, outputTimeStep);
                css.addTimeInterval(ti);
            } catch (IllegalArgumentException ex) {
                css.addTimeInterval(TimeInterval.getDefaultTimeInterval());
            }
        } else {
            for (Element e : timeIntervalElementList) {
                String s = e.getAttributeValue(XMLTags.StartTimeAttrTag);
                double startTime = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.EndTimeAttrTag);
                double endTime = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.TimeStepAttrTag);
                double timeStep = Double.valueOf(s);
                s = e.getAttributeValue(XMLTags.OutputTimeStepAttrTag);
                double outputTimeStep = Double.valueOf(s);
                TimeInterval ti = new TimeInterval(startTime, endTime, timeStep, outputTimeStep);
                css.addTimeInterval(ti);
            }
        }
        Element meshRefineElement = element.getChild(XMLTags.MeshRefinementTag, vcNamespace);
        if (meshRefineElement != null) {
            if (meshRefineElement.getChildren().size() != 0) {
                // in old model, if there is no refinement, set view level to finest
                // only set viewLevel when meshRefinement has children
                Element viewLevelChild = element.getChild(XMLTags.ViewLevelTag, vcNamespace);
                if (viewLevelChild != null) {
                    viewLevel = parseIntWithDefault(element, XMLTags.ViewLevelTag, 0);
                    css.setViewLevel(viewLevel);
                }
            }
            List<Element> levelElementList = meshRefineElement.getChildren(XMLTags.RefinementRoiTag, vcNamespace);
            for (Element levelElement : levelElementList) {
                String levelStr = levelElement.getAttributeValue(XMLTags.RefineRoiLevelAttrTag);
                int level = 1;
                if (levelStr != null) {
                    level = Integer.parseInt(levelStr);
                }
                String type = levelElement.getAttributeValue(XMLTags.RefinementRoiTypeAttrTag);
                RoiType roiType = RoiType.Membrane;
                if (type != null) {
                    try {
                        roiType = RoiType.valueOf(type);
                    } catch (Exception ex) {
                    // ignore
                    }
                }
                Element expElement = levelElement.getChild(XMLTags.ROIExpressionTag, vcNamespace);
                String roiExp = null;
                if (expElement != null) {
                    roiExp = expElement.getText();
                    RefinementRoi roi = new RefinementRoi(roiType, level, roiExp);
                    css.addRefinementRoi(roi);
                }
            }
        }
        return css;
    } catch (ExpressionException e) {
        throw new XmlParseException(e);
    }
}

Example 3 with TimeInterval

use of org.vcell.chombo.TimeInterval in project vcell by virtualcell.

the class Xmlproducer method getXML.

private Element getXML(ChomboSolverSpec chomboSolverSpec) {
    Element chomboElement = new Element(XMLTags.ChomboSolverSpec);
    Element maxBoxSize = new Element(XMLTags.MaxBoxSizeTag);
    maxBoxSize.setText(chomboSolverSpec.getMaxBoxSize() + "");
    chomboElement.addContent(maxBoxSize);
    if (!chomboSolverSpec.isFinestViewLevel()) {
        Element viewLevel = new Element(XMLTags.ViewLevelTag);
        viewLevel.setText(chomboSolverSpec.getViewLevel() + "");
        chomboElement.addContent(viewLevel);
    }
    Element fillRatio = new Element(XMLTags.FillRatioTag);
    fillRatio.setText(chomboSolverSpec.getFillRatio() + "");
    chomboElement.addContent(fillRatio);
    List<Integer> ratios = chomboSolverSpec.getRefineRatioList();
    if (ratios.size() > 0) {
        Element ratiosElement = new Element(XMLTags.RefineRatios);
        StringBuilder sb = new StringBuilder();
        for (Integer i : ratios) {
            sb.append("," + i);
        }
        ratiosElement.setText(sb.substring(1) + "");
        chomboElement.addContent(ratiosElement);
    }
    Element saveVCellOutput = new Element(XMLTags.SaveVCellOutput);
    saveVCellOutput.setText(chomboSolverSpec.isSaveVCellOutput() + "");
    chomboElement.addContent(saveVCellOutput);
    Element saveChomboOutput = new Element(XMLTags.SaveChomboOutput);
    saveChomboOutput.setText(chomboSolverSpec.isSaveChomboOutput() + "");
    chomboElement.addContent(saveChomboOutput);
    Element activateFeatureUnderDevelopment = new Element(XMLTags.ActivateFeatureUnderDevelopment);
    activateFeatureUnderDevelopment.setText(chomboSolverSpec.isActivateFeatureUnderDevelopment() + "");
    chomboElement.addContent(activateFeatureUnderDevelopment);
    Element smallVolfracThreshold = new Element(XMLTags.SmallVolfracThreshold);
    smallVolfracThreshold.setText(chomboSolverSpec.getSmallVolfracThreshold() + "");
    chomboElement.addContent(smallVolfracThreshold);
    Element blockFactor = new Element(XMLTags.BlockFactorTag);
    blockFactor.setText(chomboSolverSpec.getBlockFactor() + "");
    chomboElement.addContent(blockFactor);
    Element elementTagsGrow = new Element(XMLTags.TagsGrowTag);
    elementTagsGrow.setText(chomboSolverSpec.getTagsGrow() + "");
    chomboElement.addContent(elementTagsGrow);
    Element timeBoundsTag = new Element(XMLTags.TimeBoundTag);
    for (TimeInterval ti : chomboSolverSpec.getTimeIntervalList()) {
        Element e = new Element(XMLTags.TimeIntervalTag);
        e.setAttribute(XMLTags.StartTimeAttrTag, String.valueOf(ti.getStartingTime()));
        e.setAttribute(XMLTags.EndTimeAttrTag, String.valueOf(ti.getEndingTime()));
        e.setAttribute(XMLTags.TimeStepAttrTag, String.valueOf(ti.getTimeStep()));
        e.setAttribute(XMLTags.OutputTimeStepAttrTag, String.valueOf(ti.getOutputTimeStep()));
        timeBoundsTag.addContent(e);
    }
    chomboElement.addContent(timeBoundsTag);
    Element meshRefinement = new Element(XMLTags.MeshRefinementTag);
    for (RefinementRoi roi : chomboSolverSpec.getRefinementRois()) {
        Element roiElement = new Element(XMLTags.RefinementRoiTag);
        roiElement.setAttribute(XMLTags.RefineRoiLevelAttrTag, String.valueOf(roi.getLevel()));
        roiElement.setAttribute(XMLTags.RefinementRoiTypeAttrTag, roi.getType().name());
        if (roi.getRoiExpression() != null) {
            Element expElement = new Element(XMLTags.ROIExpressionTag);
            expElement.setText(roi.getRoiExpression().infix() + ";");
            roiElement.addContent(expElement);
        }
        meshRefinement.addContent(roiElement);
    }
    chomboElement.addContent(meshRefinement);
    return chomboElement;
}