use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class NetCDFWriter method writeHybridInputFile.
/**
* Write the model to a NetCDF file which serves as an input for stoch hybrid simulator.
* To write to a NetCDF file is a bit complicated. First, we have to create a NetCDF-3
* file. And then feed in the data.
* Creation date: (5/22/2007 5:36:03 PM)
*/
public void writeHybridInputFile(String[] parameterNames) throws Exception, cbit.vcell.parser.ExpressionException, IOException, MathException, InvalidRangeException {
Simulation simulation = simTask.getSimulation();
SimulationSymbolTable simSymbolTable = simTask.getSimulationJob().getSimulationSymbolTable();
if (initialize()) {
// we need to get model and control information first
NetcdfFileWriteable ncfile = NetcdfFileWriteable.createNew(filename, false);
// Model info. will be extracted from subDomain of mathDescription
java.util.Enumeration<SubDomain> e = simulation.getMathDescription().getSubDomains();
// remember we are dealing with compartmental model here. only 1 subdomain.
SubDomain subDomain = e.nextElement();
JumpProcess[] reactions = (JumpProcess[]) subDomain.getJumpProcesses().toArray(new JumpProcess[subDomain.getJumpProcesses().size()]);
// get species variable names
Variable[] variables = simSymbolTable.getVariables();
String[] speciesNames = new String[variables.length];
for (int i = 0; i < variables.length; i++) speciesNames[i] = variables[i].getName();
// the probabilities for reactions
Expression[] probs = new Expression[reactions.length];
for (int i = 0; i < reactions.length; i++) {
probs[i] = simSymbolTable.substituteFunctions(reactions[i].getProbabilityRate());
probs[i] = probs[i].flatten();
}
VarIniCondition[] varInis = (VarIniCondition[]) subDomain.getVarIniConditions().toArray(new VarIniCondition[subDomain.getVarIniConditions().size()]);
// the non-constant stoch variables
Vector<Variable> vars = new Vector<Variable>();
for (int i = 0; i < varInis.length; i++) {
if (varInis[i].getVar() instanceof StochVolVariable) {
vars.addElement(varInis[i].getVar());
}
}
// get reaction rate law types and rate constants
ReactionRateLaw[] reactionRateLaws = getReactionRateLaws(probs);
SolverTaskDescription solverTaskDescription = simulation.getSolverTaskDescription();
TimeBounds timeBounds = solverTaskDescription.getTimeBounds();
UniformOutputTimeSpec timeSpec = (UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec();
UniformOutputTimeSpec outputTimeSpec = ((UniformOutputTimeSpec) solverTaskDescription.getOutputTimeSpec());
NonspatialStochSimOptions stochOpt = solverTaskDescription.getStochOpt();
// create an empty NetCDF-3 file
// define dimensions
/* these sizes must match the buffers allocated in corresponding Fortran code -- see globalvariables.f90
in numerics Hy3S/src directory */
Dimension numTrial = ncfile.addDimension("NumTrials", (int) stochOpt.getNumOfTrials());
Dimension numSpecies = ncfile.addDimension("NumSpecies", vars.size());
Dimension numReactions = ncfile.addDimension("NumReactions", subDomain.getJumpProcesses().size());
int outPoints = ((int) ((timeBounds.getEndingTime() - timeBounds.getStartingTime()) / outputTimeSpec.getOutputTimeStep())) + 1;
Dimension numTimePoints = ncfile.addDimension("NumTimePoints", outPoints);
Dimension numModels = ncfile.addDimension("NumModels", 1);
Dimension numMaxDepList = ncfile.addDimension("NumMaxDepList", 6);
Dimension numMaxStoichList = ncfile.addDimension("NumMaxStoichList", 25);
Dimension stringLen = ncfile.addDimension("StringLen", 72);
// define variables
// jms info
ArrayList<Dimension> dims = new ArrayList<Dimension>();
dims.add(stringLen);
if (bMessaging) {
ncfile.addVariable("JMS_BROKER", DataType.CHAR, dims);
ncfile.addVariable("JMS_USER", DataType.CHAR, dims);
ncfile.addVariable("JMS_PASSWORD", DataType.CHAR, dims);
ncfile.addVariable("JMS_QUEUE", DataType.CHAR, dims);
ncfile.addVariable("JMS_TOPIC", DataType.CHAR, dims);
ncfile.addVariable("VCELL_USER", DataType.CHAR, dims);
ncfile.addVariable("SIMULATION_KEY", DataType.INT, new ArrayList<Dimension>());
ncfile.addVariable("JOB_INDEX", DataType.INT, new ArrayList<Dimension>());
}
// scalars
ncfile.addVariable("TStart", DataType.DOUBLE, new ArrayList<Dimension>());
ncfile.addVariable("TEnd", DataType.DOUBLE, new ArrayList<Dimension>());
ncfile.addVariable("SaveTime", DataType.DOUBLE, new ArrayList<Dimension>());
ncfile.addVariable("Volume", DataType.DOUBLE, new ArrayList<Dimension>());
ncfile.addVariable("CellGrowthTime", DataType.DOUBLE, new ArrayList<Dimension>());
ncfile.addVariable("CellGrowthTimeSD", DataType.DOUBLE, new ArrayList<Dimension>());
ncfile.addVariable("ExpType", DataType.INT, new ArrayList<Dimension>());
ncfile.addVariable("LastTrial", DataType.INT, new ArrayList<Dimension>());
ncfile.addVariable("LastModel", DataType.INT, new ArrayList<Dimension>());
ncfile.addVariable("MaxNumModels", DataType.INT, new ArrayList<Dimension>());
ncfile.addVariable("NumModels", DataType.INT, new ArrayList<Dimension>());
// variables with at least 1 dimension
ArrayList<Dimension> dimspecies = new ArrayList<Dimension>();
dimspecies.add(numSpecies);
ArrayList<Dimension> dimreactions = new ArrayList<Dimension>();
dimreactions.add(numReactions);
ncfile.addVariable("SpeciesSplitOnDivision", DataType.INT, dimspecies);
ncfile.addVariable("SaveSpeciesData", DataType.INT, dimspecies);
ncfile.addVariable("Reaction_Rate_Laws", DataType.INT, dimreactions);
ncfile.addVariable("Reaction_DListLen", DataType.INT, dimreactions);
ncfile.addVariable("Reaction_StoichListLen", DataType.INT, dimreactions);
ncfile.addVariable("Reaction_OptionalData", DataType.INT, dimreactions);
dims.clear();
dims.add(numReactions);
dims.add(numMaxStoichList);
ncfile.addVariable("Reaction_StoichCoeff", DataType.INT, dims);
ncfile.addVariable("Reaction_StoichSpecies", DataType.INT, dims);
dims.clear();
dims.add(numReactions);
dims.add(numMaxDepList);
ncfile.addVariable("Reaction_DepList", DataType.INT, dims);
dims.clear();
dims.add(numReactions);
dims.add(stringLen);
ncfile.addVariable("Reaction_names", DataType.CHAR, dims);
dims.clear();
dims.add(numSpecies);
dims.add(stringLen);
ncfile.addVariable("Species_names", DataType.CHAR, dims);
ncfile.addVariable("SpeciesIC", DataType.INT, dimspecies);
dims.clear();
dims.add(numReactions);
dims.add(numMaxDepList);
ncfile.addVariable("Reaction_Rate_Constants", DataType.DOUBLE, dims);
// create the file
try {
ncfile.create();
} catch (IOException ioe) {
ioe.printStackTrace(System.err);
throw new IOException("Error creating hybrid file " + filename + ": " + ioe.getMessage());
}
// write data to the NetCDF file
try {
// write jms info
if (bMessaging) {
ArrayChar.D1 jmsString = new ArrayChar.D1(stringLen.getLength());
String jmshost = PropertyLoader.getRequiredProperty(PropertyLoader.jmsHostExternal);
//
// Used for new REST HTTP messaging api (USE THIS WHEN Hyrbid Solvers are compiled).
//
// String jmsrestport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsRestPortExternal);
// String jmsurl = jmshost+":"+jmsrestport;
//
// connect to messaging using legacy AMQP protocol instead of new REST api. Needed for legacy pre-compiled solvers.
//
String jmsport = PropertyLoader.getRequiredProperty(PropertyLoader.jmsPortExternal);
String jmsurl = "failover:(tcp://" + jmshost + ":" + jmsport + ")";
jmsString.setString(jmsurl);
ncfile.write("JMS_BROKER", jmsString);
jmsString.setString(PropertyLoader.getRequiredProperty(PropertyLoader.jmsUser));
ncfile.write("JMS_USER", jmsString);
String jmsPassword = PropertyLoader.getSecretValue(PropertyLoader.jmsPasswordValue, PropertyLoader.jmsPasswordFile);
jmsString.setString(jmsPassword);
ncfile.write("JMS_PASSWORD", jmsString);
jmsString.setString(VCellQueue.WorkerEventQueue.getName());
ncfile.write("JMS_QUEUE", jmsString);
jmsString.setString(VCellTopic.ServiceControlTopic.getName());
ncfile.write("JMS_TOPIC", jmsString);
jmsString.setString(simulation.getVersion().getOwner().getName());
ncfile.write("VCELL_USER", jmsString);
ArrayInt.D0 scalarJMS = new ArrayInt.D0();
scalarJMS.set(Integer.parseInt(simulation.getVersion().getVersionKey() + ""));
ncfile.write("SIMULATION_KEY", scalarJMS);
scalarJMS.set(simTask.getSimulationJob().getJobIndex());
ncfile.write("JOB_INDEX", scalarJMS);
}
ArrayDouble.D0 scalarDouble = new ArrayDouble.D0();
// TStart, TEnd, SaveTime
if ((timeBounds.getEndingTime() > timeBounds.getStartingTime()) && (outputTimeSpec.getOutputTimeStep() > 0)) {
scalarDouble.set(timeBounds.getStartingTime());
ncfile.write("TStart", scalarDouble);
scalarDouble.set(timeBounds.getEndingTime());
ncfile.write("TEnd", scalarDouble);
scalarDouble.set(outputTimeSpec.getOutputTimeStep());
ncfile.write("SaveTime", scalarDouble);
} else {
System.err.println("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
throw new RuntimeException("Time setting error. Ending time smaller than starting time or save interval is not a positive value.");
}
// Volume
// we set volume to 1. This model file cannot support multi-compartmental sizes.
// When writting the rate constants, we must take the volume into account according to the reaction type.
scalarDouble.set(1);
ncfile.write("Volume", scalarDouble);
// CellGrowthTime, CellGrowthTimeSD,
scalarDouble.set(0);
ncfile.write("CellGrowthTime", scalarDouble);
ncfile.write("CellGrowthTimeSD", scalarDouble);
// ExpType, Last Trial, Last Model, MaxNumModels, NumModels
ArrayInt.D0 scalarInt = new ArrayInt.D0();
scalarInt.set(0);
ncfile.write("LastTrial", scalarInt);
ncfile.write("LastModel", scalarInt);
scalarInt.set(1);
ncfile.write("ExpType", scalarInt);
ncfile.write("MaxNumModels", scalarInt);
ncfile.write("NumModels", scalarInt);
// SpeciesSplitOnDivision
ArrayInt A1 = new ArrayInt.D1(numSpecies.getLength());
Index idx = A1.getIndex();
for (int i = 0; i < numSpecies.getLength(); i++) {
A1.setInt(idx.set(i), 0);
}
ncfile.write("SpeciesSplitOnDivision", new int[1], A1);
// SaveSpeciesData
ArrayInt A2 = new ArrayInt.D1(numSpecies.getLength());
idx = A2.getIndex();
for (int i = 0; i < numSpecies.getLength(); i++) {
A2.setInt(idx.set(i), 1);
}
ncfile.write("SaveSpeciesData", new int[1], A2);
// Reaction_Rate_Laws
ArrayInt A3 = new ArrayInt.D1(numReactions.getLength());
idx = A3.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
A3.setInt(idx.set(i), reactionRateLaws[i].getLawType());
}
ncfile.write("Reaction_Rate_Laws", new int[1], A3);
// Reaction_DListLen
ArrayInt A4 = new ArrayInt.D1(numReactions.getLength());
idx = A4.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_0)
A4.setInt(idx.set(i), 0);
else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_1) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_1substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_1substrate))
A4.setInt(idx.set(i), 1);
else if ((reactionRateLaws[i].getLawType() == ReactionRateLaw.order_2_2substrate) || (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_2substrate))
A4.setInt(idx.set(i), 2);
else if (reactionRateLaws[i].getLawType() == ReactionRateLaw.order_3_3substrate)
A4.setInt(idx.set(i), 3);
}
ncfile.write("Reaction_DListLen", new int[1], A4);
// Reaction_StoichListLen
ArrayInt A5 = new ArrayInt.D1(numReactions.getLength());
idx = A5.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
A5.setInt(idx.set(i), reactions[i].getActions().size());
}
ncfile.write("Reaction_StoichListLen", new int[1], A5);
// Reaction_OptionalData
ArrayInt A6 = new ArrayInt.D1(numReactions.getLength());
idx = A6.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
A6.setInt(idx.set(i), 0);
}
ncfile.write("Reaction_OptionalData", new int[1], A6);
// Reaction_StoichCoeff
ArrayInt A7 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
idx = A7.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
Action[] actions = (Action[]) reactions[i].getActions().toArray(new Action[reactions[i].getActions().size()]);
for (int j = 0; j < actions.length; j++) {
try {
actions[j].getOperand().evaluateConstant();
int coeff = (int) Math.round(actions[j].getOperand().evaluateConstant());
A7.setInt(idx.set(i, j), coeff);
} catch (ExpressionException ex) {
ex.printStackTrace(System.err);
throw new ExpressionException(ex.getMessage());
}
}
}
ncfile.write("Reaction_StoichCoeff", new int[2], A7);
// Reaction_StoichSpecies
ArrayInt A8 = new ArrayInt.D2(numReactions.getLength(), numMaxStoichList.getLength());
idx = A8.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
ArrayList<Action> actions = reactions[i].getActions();
for (int j = 0; j < actions.size(); j++) {
A8.setInt(idx.set(i, j), getVariableIndex(((Action) actions.get(j)).getVar().getName(), vars));
}
}
ncfile.write("Reaction_StoichSpecies", new int[2], A8);
// Reaction_DepList
ArrayInt A9 = new ArrayInt.D2(numReactions.getLength(), numMaxDepList.getLength());
idx = A9.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
ReactionRateLaw rl = reactionRateLaws[i];
Hashtable<String, Integer> tem = varInProbOrderHash[i];
Enumeration<String> varnames = tem.keys();
if (rl.getLawType() == ReactionRateLaw.order_0) {
// don't do anything here.
} else if ((rl.getLawType() == ReactionRateLaw.order_1) || (rl.getLawType() == ReactionRateLaw.order_2_1substrate) || (rl.getLawType() == ReactionRateLaw.order_3_1substrate) || (rl.getLawType() == ReactionRateLaw.order_2_2substrate) || (rl.getLawType() == ReactionRateLaw.order_3_3substrate)) {
int j = 0;
while (varnames.hasMoreElements()) {
String name = varnames.nextElement();
A9.setInt(idx.set(i, j), getVariableIndex(name, vars));
j++;
}
} else if (rl.getLawType() == ReactionRateLaw.order_3_2substrate) {
int order = 0;
String highOrderName = "";
String lowOrderName = "";
// we must make sure to put the higher order species first.
while (varnames.hasMoreElements()) {
lowOrderName = varnames.nextElement();
if (tem.get(lowOrderName) > order) {
String s = highOrderName;
highOrderName = lowOrderName;
lowOrderName = s;
order = tem.get(highOrderName);
}
}
A9.setInt(idx.set(i, 0), getVariableIndex(highOrderName, vars));
A9.setInt(idx.set(i, 1), getVariableIndex(lowOrderName, vars));
}
}
ncfile.write("Reaction_DepList", new int[2], A9);
// Reaction_names
ArrayChar A10 = new ArrayChar.D2(numReactions.getLength(), stringLen.getLength());
for (int i = 0; i < numReactions.getLength(); i++) {
String name = reactions[i].getName();
int diff = stringLen.getLength() - name.length();
if (diff >= 0) {
for (int j = 0; j < diff; j++) {
name = name + " ";
}
A10.setString(i, name);
} else
throw new RuntimeException("Name of Reaction:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
}
ncfile.write("Reaction_names", A10);
// Species_names
ArrayChar A11 = new ArrayChar.D2(numSpecies.getLength(), stringLen.getLength());
for (int i = 0; i < numSpecies.getLength(); i++) {
String name = vars.elementAt(i).getName();
int diff = stringLen.getLength() - name.length();
if (diff >= 0) {
for (int j = 0; j < diff; j++) {
name = name + " ";
}
A11.setString(i, name);
} else
throw new RuntimeException("Name of Species:" + name + " is too long. Please shorten to " + stringLen.getLength() + " chars.");
}
ncfile.write("Species_names", A11);
// Species Initial Condition (in number of molecules).
// Species iniCondition are sampled from a poisson distribution(which has a mean of the current iniExp value)
RandomDataGenerator dist = new RandomDataGenerator();
if (stochOpt.isUseCustomSeed()) {
Integer randomSeed = stochOpt.getCustomSeed();
if (randomSeed != null) {
dist.reSeed(randomSeed);
}
}
ArrayLong A12 = new ArrayLong.D1(numSpecies.getLength());
idx = A12.getIndex();
for (int i = 0; i < numSpecies.getLength(); i++) {
try {
VarIniCondition varIniCondition = subDomain.getVarIniCondition(vars.elementAt(i));
Expression varIniExp = varIniCondition.getIniVal();
varIniExp.bindExpression(simSymbolTable);
varIniExp = simSymbolTable.substituteFunctions(varIniExp).flatten();
double expectedCount = varIniExp.evaluateConstant();
long varCount = 0;
if (varIniCondition instanceof VarIniCount) {
varCount = (long) expectedCount;
} else {
if (expectedCount > 0) {
varCount = dist.nextPoisson(expectedCount);
}
}
A12.setLong(idx.set(i), varCount);
} catch (ExpressionException ex) {
ex.printStackTrace(System.err);
throw new ExpressionException(ex.getMessage());
}
}
ncfile.write("SpeciesIC", new int[1], A12);
// Reaction_Rate_Constants(NumReactions, NumMaxDepList) ;
ArrayDouble A13 = new ArrayDouble.D2(numReactions.getLength(), numMaxDepList.getLength());
idx = A13.getIndex();
for (int i = 0; i < numReactions.getLength(); i++) {
ReactionRateLaw rl = reactionRateLaws[i];
A13.setDouble(idx.set(i, 0), rl.getRateConstant());
}
ncfile.write("Reaction_Rate_Constants", A13);
} catch (IOException ioe) {
ioe.printStackTrace(System.err);
throw new IOException("Error writing hybrid input file " + filename + ": " + ioe.getMessage());
} catch (InvalidRangeException ire) {
ire.printStackTrace(System.err);
throw new InvalidRangeException("Error writing hybrid input file " + filename + ": " + ire.getMessage());
}
try {
ncfile.close();
} catch (IOException ioe) {
throw new IOException("Error closing file " + filename + ". " + ioe.getMessage());
}
}
}
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class RunSims method isSmoldynTimeStepOK.
private boolean isSmoldynTimeStepOK(Simulation sim) {
for (int jobIndex = 0; jobIndex < sim.getScanCount(); jobIndex++) {
SimulationSymbolTable simSymbolTable = new SimulationSymbolTable(sim, jobIndex);
double Dmax = 0;
MathDescription mathDesc = sim.getMathDescription();
Enumeration<SubDomain> subDomainEnumeration = mathDesc.getSubDomains();
while (subDomainEnumeration.hasMoreElements()) {
SubDomain subDomain = subDomainEnumeration.nextElement();
// }
for (ParticleProperties particleProperties : subDomain.getParticleProperties()) {
try {
Expression newExp = new Expression(particleProperties.getDiffusion());
newExp.bindExpression(simSymbolTable);
newExp = simSymbolTable.substituteFunctions(newExp).flatten();
try {
double diffConstant = newExp.evaluateConstant();
Dmax = Math.max(Dmax, diffConstant);
} catch (ExpressionException ex) {
throw new ExpressionException("diffusion coefficient for variable " + particleProperties.getVariable().getQualifiedName() + " is not a constant. Constants are required for all diffusion coefficients");
}
} catch (Exception ex) {
}
}
}
double s = sim.getMeshSpecification().getDx(sim.hasCellCenteredMesh());
double dt = sim.getSolverTaskDescription().getTimeStep().getDefaultTimeStep();
if (dt >= s * s / (2 * Dmax)) {
smoldynTimestepVars = new SmoldynTimeStepVars(s, Dmax);
return false;
}
}
return true;
}
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class SimulationData method getChomboFiles.
@Override
public ChomboFiles getChomboFiles() throws IOException, XmlParseException, ExpressionException {
if (chomboFileIterationIndices == null) {
throw new RuntimeException("SimulationData.chomboFileIterationIndices is null, can't process Chombo HDF5 files");
}
if (!(getVcDataId() instanceof VCSimulationDataIdentifier)) {
throw new RuntimeException("SimulationData.getVcDataId() is not a VCSimulationDataIdentifier (type is " + getVcDataId().getClass().getName() + "), can't process chombo HDF5 files");
}
VCSimulationDataIdentifier vcDataID = (VCSimulationDataIdentifier) getVcDataId();
String expectedMeshfile = vcDataID.getID() + ".mesh.hdf5";
File meshFile = amplistorHelper.getFile(expectedMeshfile);
ChomboFiles chomboFiles = new ChomboFiles(vcDataID.getSimulationKey(), vcDataID.getJobIndex(), meshFile);
String simtaskFilePath = vcDataID.getID() + "_0.simtask.xml";
File simtaskFile = amplistorHelper.getFile(simtaskFilePath);
if (!simtaskFile.exists()) {
throw new RuntimeException("Chombo dataset mission .simtask.xml file, please rerun");
}
String xmlString = FileUtils.readFileToString(simtaskFile);
SimulationTask simTask = XmlHelper.XMLToSimTask(xmlString);
if (!simTask.getSimulation().getSolverTaskDescription().getChomboSolverSpec().isSaveChomboOutput() && !simTask.getSimulation().getSolverTaskDescription().isParallel()) {
throw new RuntimeException("Export of Chombo simulations to VTK requires chombo data, select 'Chombo' data format in simulation solver options and rerun simulation.");
}
CartesianMeshChombo chomboMesh = (CartesianMeshChombo) mesh;
FeaturePhaseVol[] featurePhaseVols = chomboMesh.getFeaturePhaseVols();
for (int timeIndex : chomboFileIterationIndices) {
if (featurePhaseVols == null) {
// for old format which doesn't have featurephasevols, we need to try ivol up to 20 for each feature
Enumeration<SubDomain> subdomainEnum = simTask.getSimulation().getMathDescription().getSubDomains();
while (subdomainEnum.hasMoreElements()) {
SubDomain subDomain = subdomainEnum.nextElement();
if (subDomain instanceof CompartmentSubDomain) {
for (int ivol = 0; ivol < 20; ++ivol) {
// can be many vol, let us try 20
findChomboFeatureVolFile(chomboFiles, vcDataID, subDomain.getName(), ivol, timeIndex);
}
}
}
} else {
// note: some feature + ivol doesn't have a file if there are no variables defined in that feature
for (FeaturePhaseVol pfv : featurePhaseVols) {
findChomboFeatureVolFile(chomboFiles, vcDataID, pfv.feature, pfv.ivol, timeIndex);
}
}
}
return chomboFiles;
}
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class ModelOptimizationSpec method fromMath.
/**
* Insert the method's description here.
* Creation date: (11/1/2005 8:11:32 PM)
* @return cbit.vcell.mapping.MathSystemHash
* @param mathDesc cbit.vcell.math.MathDescription
*/
private static MathSystemHash fromMath(cbit.vcell.math.MathDescription mathDesc) {
MathSystemHash hash = new MathSystemHash();
hash.addSymbol(new MathSystemHash.IndependentVariable("t"));
hash.addSymbol(new MathSystemHash.IndependentVariable("x"));
hash.addSymbol(new MathSystemHash.IndependentVariable("y"));
hash.addSymbol(new MathSystemHash.IndependentVariable("z"));
Variable[] vars = (Variable[]) BeanUtils.getArray(mathDesc.getVariables(), Variable.class);
for (int i = 0; i < vars.length; i++) {
hash.addSymbol(new MathSystemHash.Variable(vars[i].getName(), vars[i].getExpression()));
}
SubDomain[] subDomains = (SubDomain[]) BeanUtils.getArray(mathDesc.getSubDomains(), SubDomain.class);
for (int i = 0; i < subDomains.length; i++) {
Equation[] equations = (Equation[]) BeanUtils.getArray(subDomains[i].getEquations(), Equation.class);
for (int j = 0; j < equations.length; j++) {
MathSystemHash.Variable var = (MathSystemHash.Variable) hash.getSymbol(equations[j].getVariable().getName());
hash.addSymbol(new MathSystemHash.VariableInitial(var, equations[j].getInitialExpression()));
if (equations[j] instanceof PdeEquation) {
// hash.addSymbol(new MathSystemHash.VariableDerivative(var,pde.getRateExpression()));
throw new RuntimeException("MathSystemHash doesn't yet support spatial models");
} else if (equations[j] instanceof VolumeRegionEquation) {
// hash.addSymbol(new MathSystemHash.VariableDerivative(var,vre.getRateExpression()));
throw new RuntimeException("MathSystemHash doesn't yet support spatial models");
} else if (equations[j] instanceof MembraneRegionEquation) {
// hash.addSymbol(new MathSystemHash.VariableDerivative(var,mre.getRateExpression()));
throw new RuntimeException("MathSystemHash doesn't yet support spatial models");
} else if (equations[j] instanceof FilamentRegionEquation) {
// hash.addSymbol(new MathSystemHash.VariableDerivative(var,fre.getRateExpression()));
throw new RuntimeException("MathSystemHash doesn't yet support spatial models");
} else if (equations[j] instanceof OdeEquation) {
OdeEquation ode = (OdeEquation) equations[j];
hash.addSymbol(new MathSystemHash.VariableDerivative(var, ode.getRateExpression()));
}
}
}
return hash;
}
use of cbit.vcell.math.SubDomain in project vcell by virtualcell.
the class ComsolModelBuilder method getVCCModel.
public static VCCModel getVCCModel(SimulationJob vcellSimJob) throws ExpressionException {
MathDescription vcellMathDesc = vcellSimJob.getSimulation().getMathDescription();
Geometry vcellGeometry = vcellMathDesc.getGeometry();
GeometrySpec vcellGeometrySpec = vcellGeometry.getGeometrySpec();
int vcellDim = vcellGeometrySpec.getDimension();
VCCModel model = new VCCModel("Model", vcellDim);
model.modelpath = "D:\\Developer\\eclipse\\workspace_refactor\\comsol_java\\src";
model.comments = "Untitled\n\n";
VCCModelNode comp1 = new VCCModelNode("comp1");
model.modelnodes.add(comp1);
// if (vcellDim != 2){
// throw new RuntimeException("expecting 2D simulation");
// }
//
// assume initial geometry is circle centered at 0.5, 0.5 of radius 0.3
//
// String comsolOutsideDomainName = "dif1";
// String comsolInsideDomainName = "c1";
VCCGeomSequence geom1 = new VCCGeomSequence("geom1", vcellDim);
model.geometrysequences.add(geom1);
VCCMeshSequence mesh1 = new VCCMeshSequence("mesh1", geom1);
model.meshes.add(mesh1);
VCCStudy std1 = new VCCStudy("std1");
model.study = std1;
TimeBounds timeBounds = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeBounds();
TimeStep timeStep = vcellSimJob.getSimulation().getSolverTaskDescription().getTimeStep();
String beginTime = Double.toString(timeBounds.getStartingTime());
String endTime = Double.toString(timeBounds.getEndingTime());
String step = Double.toString(timeStep.getDefaultTimeStep());
VCCStudyFeature time = new VCCTransientStudyFeature("time", beginTime, step, endTime);
std1.features.add(time);
if (vcellGeometrySpec.getImage() != null) {
throw new RuntimeException("image-based geometries not yet supported by VCell's COMSOL model builder");
}
if (vcellGeometrySpec.getNumSubVolumes() == 0) {
throw new RuntimeException("no subvolumes defined in geometry");
}
if (vcellGeometrySpec.getNumAnalyticOrCSGSubVolumes() != vcellGeometrySpec.getNumSubVolumes()) {
throw new RuntimeException("only analytic and CSG subvolumes currently supported by VCell's COMSOL model builder");
}
//
// add geometry for all subvolumes
//
HashMap<String, VCCGeomFeature> subvolumeNameFeatureMap = new HashMap<String, VCCGeomFeature>();
SubVolume[] subVolumes = vcellGeometrySpec.getSubVolumes();
for (int i = 0; i < subVolumes.length; i++) {
SubVolume subvolume = subVolumes[i];
if (subvolume instanceof CSGObject) {
CSGObject vcellCSGObject = (CSGObject) subvolume;
CSGNode vcellCSGNode = vcellCSGObject.getRoot();
ArrayList<VCCGeomFeature> geomFeatureList = new ArrayList<VCCGeomFeature>();
VCCGeomFeature feature = csgVisitor(vcellCSGNode, geomFeatureList, subvolume.getName());
geom1.geomfeatures.addAll(geomFeatureList);
if (i == 0) {
// first subvolume (on top in ordinals) doesn't need any differencing
subvolumeNameFeatureMap.put(subvolume.getName(), feature);
} else {
// have to subtract union of prior subvolumes
ArrayList<VCCGeomFeature> priorFeatures = new ArrayList<VCCGeomFeature>();
for (int j = 0; j < i; j++) {
CSGObject priorCSGObject = (CSGObject) subVolumes[j];
CSGNode priorCSGNode = priorCSGObject.getRoot();
geomFeatureList.clear();
VCCGeomFeature priorFeature = csgVisitor(priorCSGNode, geomFeatureList, subvolume.getName());
priorFeatures.add(priorFeature);
geom1.geomfeatures.addAll(geomFeatureList);
}
VCCDifference diff = new VCCDifference("diff" + subvolume.getName(), Keep.off);
diff.input.add(feature);
diff.input2.addAll(priorFeatures);
geom1.geomfeatures.add(diff);
subvolumeNameFeatureMap.put(subvolume.getName(), diff);
}
} else {
throw new RuntimeException("only CSG subvolumes currently supported by VCell's COMSOL model builder");
}
}
//
// add geometry for all surfaceClasses
//
HashMap<String, VCCGeomFeature> surfaceclassNameFeatureMap = new HashMap<String, VCCGeomFeature>();
SurfaceClass[] surfaceClasses = vcellGeometry.getGeometrySurfaceDescription().getSurfaceClasses();
for (int i = 0; i < surfaceClasses.length; i++) {
SurfaceClass surfaceClass = surfaceClasses[i];
Set<SubVolume> adjacentSubvolumes = surfaceClass.getAdjacentSubvolumes();
if (adjacentSubvolumes.size() != 2) {
throw new RuntimeException("expecting two adjacent subvolumes for surface " + surfaceClass.getName() + " in COMSOL model builder");
}
// find adjacent Geometry Features (for subvolumes)
Iterator<SubVolume> svIter = adjacentSubvolumes.iterator();
SubVolume subvolume0 = svIter.next();
SubVolume subvolume1 = svIter.next();
ArrayList<VCCGeomFeature> adjacentFeatures = new ArrayList<VCCGeomFeature>();
adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume0.getName()));
adjacentFeatures.add(subvolumeNameFeatureMap.get(subvolume1.getName()));
String name = "inter_" + subvolume0.getName() + "_" + subvolume1.getName();
// surfaces are dimension N-1
int entitydim = vcellDim - 1;
VCCIntersectionSelection intersect_subvolumes = new VCCIntersectionSelection(name, entitydim);
intersect_subvolumes.input.addAll(adjacentFeatures);
geom1.geomfeatures.add(intersect_subvolumes);
surfaceclassNameFeatureMap.put(surfaceClass.getName(), intersect_subvolumes);
}
SimulationSymbolTable symbolTable = new SimulationSymbolTable(vcellSimJob.getSimulation(), vcellSimJob.getJobIndex());
//
for (SubDomain subDomain : Collections.list(vcellMathDesc.getSubDomains())) {
for (Equation equ : subDomain.getEquationCollection()) {
if (equ instanceof PdeEquation || equ instanceof OdeEquation) {
VCCGeomFeature geomFeature = null;
final int dim;
if (subDomain instanceof CompartmentSubDomain) {
geomFeature = subvolumeNameFeatureMap.get(subDomain.getName());
dim = vcellDim;
} else if (subDomain instanceof MembraneSubDomain) {
geomFeature = surfaceclassNameFeatureMap.get(subDomain.getName());
dim = vcellDim - 1;
} else {
throw new RuntimeException("subdomains of type '" + subDomain.getClass().getSimpleName() + "' not yet supported in COMSOL model builder");
}
if (geomFeature == null) {
throw new RuntimeException("cannot find COMSOL geometry feature named " + subDomain.getName() + " in COMSOL model builder");
}
VCCConvectionDiffusionEquation cdeq = new VCCConvectionDiffusionEquation("cdeq_" + equ.getVariable().getName(), geom1, geomFeature, dim);
cdeq.fieldName = equ.getVariable().getName();
cdeq.initial = MathUtilities.substituteModelParameters(equ.getInitialExpression(), symbolTable).flatten().infix();
cdeq.sourceTerm_f = MathUtilities.substituteModelParameters(equ.getRateExpression(), symbolTable).flatten().infix();
if (equ instanceof PdeEquation) {
PdeEquation pde = (PdeEquation) equ;
cdeq.diffTerm_c = MathUtilities.substituteModelParameters(pde.getDiffusionExpression(), symbolTable).flatten().infix();
if (subDomain instanceof CompartmentSubDomain) {
CompartmentSubDomain compartmentSubdomain = (CompartmentSubDomain) subDomain;
ArrayList<String> be = new ArrayList<String>();
if (pde.getVelocityX() != null) {
be.add(MathUtilities.substituteModelParameters(pde.getVelocityX(), symbolTable).flatten().infix());
} else {
be.add("0");
}
if (vcellDim >= 2) {
if (pde.getVelocityY() != null) {
be.add(MathUtilities.substituteModelParameters(pde.getVelocityY(), symbolTable).flatten().infix());
} else {
be.add("0");
}
}
if (vcellDim == 3) {
if (pde.getVelocityY() != null) {
be.add(MathUtilities.substituteModelParameters(pde.getVelocityZ(), symbolTable).flatten().infix());
} else {
be.add("0");
}
}
cdeq.advection_be = be.toArray(new String[vcellDim]);
//
// look for membrane boundary conditions for this variable
//
MembraneSubDomain[] membraneSubdomains = vcellMathDesc.getMembraneSubDomains(compartmentSubdomain);
for (MembraneSubDomain membraneSubdomain : membraneSubdomains) {
JumpCondition jumpCondition = membraneSubdomain.getJumpCondition((VolVariable) pde.getVariable());
if (jumpCondition != null) {
Expression fluxExpr = null;
if (membraneSubdomain.getInsideCompartment() == compartmentSubdomain) {
fluxExpr = jumpCondition.getInFluxExpression();
} else if (membraneSubdomain.getOutsideCompartment() == compartmentSubdomain) {
fluxExpr = jumpCondition.getOutFluxExpression();
}
String name = equ.getVariable().getName() + "_flux_" + membraneSubdomain.getName();
VCCGeomFeature selection = surfaceclassNameFeatureMap.get(membraneSubdomain.getName());
VCCFluxBoundary fluxBoundary = new VCCFluxBoundary(name, selection, vcellDim - 1);
fluxBoundary.flux_g = MathUtilities.substituteModelParameters(fluxExpr, symbolTable).flatten().infix();
cdeq.features.add(fluxBoundary);
}
}
}
}
model.physics.add(cdeq);
}
}
}
//
return model;
}
Aggregations