Search in sources :

Example 26 with MembraneSubDomain

use of cbit.vcell.math.MembraneSubDomain 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;
}
Also used : JumpCondition(cbit.vcell.math.JumpCondition) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) MathDescription(cbit.vcell.math.MathDescription) HashMap(java.util.HashMap) SurfaceClass(cbit.vcell.geometry.SurfaceClass) VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation) CSGNode(cbit.vcell.geometry.CSGNode) ArrayList(java.util.ArrayList) VCCFluxBoundary(org.vcell.solver.comsol.model.VCCPhysicsFeature.VCCFluxBoundary) GeometrySpec(cbit.vcell.geometry.GeometrySpec) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) PdeEquation(cbit.vcell.math.PdeEquation) TimeBounds(cbit.vcell.solver.TimeBounds) TimeStep(cbit.vcell.solver.TimeStep) SubVolume(cbit.vcell.geometry.SubVolume) VCCStudyFeature(org.vcell.solver.comsol.model.VCCStudyFeature) CSGObject(cbit.vcell.geometry.CSGObject) VCCGeomFeature(org.vcell.solver.comsol.model.VCCGeomFeature) VCCDifference(org.vcell.solver.comsol.model.VCCGeomFeature.VCCDifference) VCCGeomSequence(org.vcell.solver.comsol.model.VCCGeomSequence) VCCIntersectionSelection(org.vcell.solver.comsol.model.VCCGeomFeature.VCCIntersectionSelection) SimulationSymbolTable(cbit.vcell.solver.SimulationSymbolTable) VCCConvectionDiffusionEquation(org.vcell.solver.comsol.model.VCCConvectionDiffusionEquation) OdeEquation(cbit.vcell.math.OdeEquation) PdeEquation(cbit.vcell.math.PdeEquation) Equation(cbit.vcell.math.Equation) VCCModelNode(org.vcell.solver.comsol.model.VCCModelNode) VCCTransientStudyFeature(org.vcell.solver.comsol.model.VCCTransientStudyFeature) VCCMeshSequence(org.vcell.solver.comsol.model.VCCMeshSequence) Geometry(cbit.vcell.geometry.Geometry) VCCStudy(org.vcell.solver.comsol.model.VCCStudy) OdeEquation(cbit.vcell.math.OdeEquation) Expression(cbit.vcell.parser.Expression) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) VCCModel(org.vcell.solver.comsol.model.VCCModel)

Example 27 with MembraneSubDomain

use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.

the class ParticleMathMapping method combineHybrid.

private void combineHybrid() throws MappingException, ExpressionException, MatrixException, MathException, ModelException {
    ArrayList<SpeciesContext> continuousSpecies = new ArrayList<SpeciesContext>();
    ArrayList<ParticleVariable> continuousSpeciesParticleVars = new ArrayList<ParticleVariable>();
    ArrayList<SpeciesContext> stochSpecies = new ArrayList<SpeciesContext>();
    // 
    // categorize speciesContexts as continuous and stochastic
    // 
    SpeciesContextSpec[] scsArray = getSimulationContext().getReactionContext().getSpeciesContextSpecs();
    continuousSpecies = new ArrayList<SpeciesContext>();
    stochSpecies = new ArrayList<SpeciesContext>();
    for (SpeciesContextSpec speciesContextSpec : scsArray) {
        if (!getSimulationContext().isStoch() || speciesContextSpec.isForceContinuous()) {
            continuousSpecies.add(speciesContextSpec.getSpeciesContext());
            Variable variable = getMathSymbolMapping().getVariable(speciesContextSpec.getSpeciesContext());
            if (variable instanceof ParticleVariable) {
                continuousSpeciesParticleVars.add((ParticleVariable) variable);
            }
        } else {
            stochSpecies.add(speciesContextSpec.getSpeciesContext());
        }
    }
    if (continuousSpecies.isEmpty()) {
        return;
    }
    // 
    // create continuous mathDescription ... add stochastic variables and processes to the continuous Math and use this.
    // 
    DiffEquMathMapping mathMapping = new DiffEquMathMapping(getSimulationContext(), callback, networkGenerationRequirements);
    mathMapping.refresh(null);
    MathDescription contMathDesc = mathMapping.getMathDescription();
    // 
    // get list of all continuous variables
    // 
    HashMap<String, Variable> allContinuousVars = new HashMap<String, Variable>();
    Enumeration<Variable> enumVar = contMathDesc.getVariables();
    while (enumVar.hasMoreElements()) {
        Variable var = enumVar.nextElement();
        allContinuousVars.put(var.getName(), var);
    }
    // 
    // replace those continuous variables and equations for stochastic speciesContexts
    // with the particleVariables and particleProperties
    // (ParticleJumpProcesses removed later)
    // 
    ModelUnitSystem unitSystem = getSimulationContext().getModel().getUnitSystem();
    for (SpeciesContext stochSpeciesContext : stochSpecies) {
        Variable contVar = mathMapping.getMathSymbolMapping().getVariable(stochSpeciesContext);
        Variable stochVar = getMathSymbolMapping().getVariable(stochSpeciesContext);
        allContinuousVars.put(stochVar.getName(), stochVar);
        // 
        // replace continuous "concentration" VolVariable/MemVariable for this particle with a Function for concentration
        // 
        allContinuousVars.remove(contVar);
        VCUnitDefinition sizeUnit = unitSystem.getLengthUnit().raiseTo(new RationalNumber(stochSpeciesContext.getStructure().getDimension()));
        VCUnitDefinition stochasticDensityUnit = unitSystem.getStochasticSubstanceUnit().divideBy(sizeUnit);
        VCUnitDefinition continuousDensityUnit = unitSystem.getConcentrationUnit(stochSpeciesContext.getStructure());
        if (stochasticDensityUnit.isEquivalent(continuousDensityUnit)) {
            allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), new Expression(stochVar, getNameScope()), contVar.getDomain()));
        } else {
            Expression conversionFactorExp = getUnitFactor(continuousDensityUnit.divideBy(stochasticDensityUnit));
            allContinuousVars.put(contVar.getName(), new Function(contVar.getName(), Expression.mult(new Expression(stochVar, getNameScope()), conversionFactorExp), contVar.getDomain()));
        }
        // 
        // remove continuous equation
        // 
        Enumeration<SubDomain> contSubDomains = contMathDesc.getSubDomains();
        while (contSubDomains.hasMoreElements()) {
            SubDomain contSubDomain = contSubDomains.nextElement();
            contSubDomain.removeEquation(contVar);
            if (contSubDomain instanceof MembraneSubDomain) {
                ((MembraneSubDomain) contSubDomain).removeJumpCondition(contVar);
            }
        }
        // 
        // remove all continuous variables for speciesContextSpec parameters (e.g. initial conditions, diffusion rates, boundary conditions, velocities)
        // 
        SpeciesContextSpec scs = getSimulationContext().getReactionContext().getSpeciesContextSpec(stochSpeciesContext);
        Parameter[] scsParameters = scs.getParameters();
        for (Parameter parameter : scsParameters) {
            Variable continuousScsParmVariable = mathMapping.getMathSymbolMapping().getVariable(parameter);
            allContinuousVars.remove(continuousScsParmVariable);
        }
        // 
        // copy ParticleJumpProcess and ParticleProperties to the continuous math
        // 
        SubDomain contSubDomain = contMathDesc.getSubDomain(contVar.getDomain().getName());
        SubDomain stochSubDomain = mathDesc.getSubDomain(stochVar.getDomain().getName());
        ParticleProperties particleProperties = stochSubDomain.getParticleProperties(stochVar);
        contSubDomain.addParticleProperties(particleProperties);
    }
    // 
    // add all ParticleJumpProcesses to the continuous model
    // 
    Enumeration<SubDomain> enumStochSubdomains = mathDesc.getSubDomains();
    while (enumStochSubdomains.hasMoreElements()) {
        SubDomain stochSubdomain = enumStochSubdomains.nextElement();
        SubDomain contSubdomain = contMathDesc.getSubDomain(stochSubdomain.getName());
        for (ParticleJumpProcess particleJumpProcess : stochSubdomain.getParticleJumpProcesses()) {
            // 
            // modify "selection list" (particleVariables), probability rate, and actions if referenced particleVariable is to be "forced continuous"
            // 
            ParticleVariable[] selectedParticles = particleJumpProcess.getParticleVariables();
            for (ParticleVariable particleVariable : selectedParticles) {
                if (continuousSpeciesParticleVars.contains(particleVariable)) {
                    particleJumpProcess.remove(particleVariable);
                    JumpProcessRateDefinition jumpProcessRateDefinition = particleJumpProcess.getParticleRateDefinition();
                    if (jumpProcessRateDefinition instanceof MacroscopicRateConstant) {
                        MacroscopicRateConstant macroscopicRateConstant = (MacroscopicRateConstant) jumpProcessRateDefinition;
                        macroscopicRateConstant.setExpression(Expression.mult(macroscopicRateConstant.getExpression(), new Expression(particleVariable, null)));
                    } else if (jumpProcessRateDefinition instanceof InteractionRadius) {
                        throw new MappingException("cannot adjust interaction radius for reaction process " + particleJumpProcess.getName() + ", particle " + particleVariable.getName() + " is continuous");
                    } else {
                        throw new MappingException("rate definition type " + jumpProcessRateDefinition.getClass().getSimpleName() + " not yet implemented for hybrid PDE/Particle math generation");
                    }
                }
                Iterator<Action> iterAction = particleJumpProcess.getActions().iterator();
                while (iterAction.hasNext()) {
                    Action action = iterAction.next();
                    if (continuousSpeciesParticleVars.contains(action.getVar())) {
                        iterAction.remove();
                    }
                }
            }
            if (!particleJumpProcess.getActions().isEmpty()) {
                contSubdomain.addParticleJumpProcess(particleJumpProcess);
            }
        }
    }
    // 
    for (MathMappingParameter mathMappingParameter : fieldMathMappingParameters) {
        if (mathMappingParameter instanceof UnitFactorParameter) {
            String name = mathMappingParameter.getName();
            if (!allContinuousVars.containsKey(name)) {
                allContinuousVars.put(name, newFunctionOrConstant(name, mathMappingParameter.getExpression(), null));
            }
        }
    }
    // 
    // add constants and functions from the particle math that aren't already defined in the continuous math
    // 
    Enumeration<Variable> enumVars = mathDesc.getVariables();
    while (enumVars.hasMoreElements()) {
        Variable var = enumVars.nextElement();
        if (var instanceof Constant || var instanceof Function) {
            String name = var.getName();
            if (!allContinuousVars.containsKey(name)) {
                allContinuousVars.put(name, var);
            }
        }
    }
    contMathDesc.setAllVariables(allContinuousVars.values().toArray(new Variable[0]));
    mathDesc = contMathDesc;
    // 
    for (int i = 0; i < fieldMathMappingParameters.length; i++) {
        if (fieldMathMappingParameters[i] instanceof UnitFactorParameter) {
            GeometryClass geometryClass = fieldMathMappingParameters[i].getGeometryClass();
            Variable variable = newFunctionOrConstant(getMathSymbol(fieldMathMappingParameters[i], geometryClass), getIdentifierSubstitutions(fieldMathMappingParameters[i].getExpression(), fieldMathMappingParameters[i].getUnitDefinition(), geometryClass), fieldMathMappingParameters[i].getGeometryClass());
            if (mathDesc.getVariable(variable.getName()) == null) {
                mathDesc.addVariable(variable);
            }
        }
    }
    if (!mathDesc.isValid()) {
        System.out.println(mathDesc.getVCML_database());
        throw new MappingException("generated an invalid mathDescription: " + mathDesc.getWarning());
    }
    System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string begin ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
    System.out.println(mathDesc.getVCML());
    System.out.println("]]]]]]]]]]]]]]]]]]]]]] VCML string end ]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]");
}
Also used : GeometryClass(cbit.vcell.geometry.GeometryClass) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Action(cbit.vcell.math.Action) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable) ParticleVariable(cbit.vcell.math.ParticleVariable) Variable(cbit.vcell.math.Variable) MathDescription(cbit.vcell.math.MathDescription) HashMap(java.util.HashMap) VolumeParticleVariable(cbit.vcell.math.VolumeParticleVariable) MembraneParticleVariable(cbit.vcell.math.MembraneParticleVariable) ParticleVariable(cbit.vcell.math.ParticleVariable) MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant) Constant(cbit.vcell.math.Constant) ArrayList(java.util.ArrayList) SpeciesContext(cbit.vcell.model.SpeciesContext) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Function(cbit.vcell.math.Function) MacroscopicRateConstant(cbit.vcell.math.MacroscopicRateConstant) RationalNumber(ucar.units_vcell.RationalNumber) ModelUnitSystem(cbit.vcell.model.ModelUnitSystem) JumpProcessRateDefinition(cbit.vcell.math.JumpProcessRateDefinition) InteractionRadius(cbit.vcell.math.InteractionRadius) ParticleJumpProcess(cbit.vcell.math.ParticleJumpProcess) VCUnitDefinition(cbit.vcell.units.VCUnitDefinition) Expression(cbit.vcell.parser.Expression) StructureMappingParameter(cbit.vcell.mapping.StructureMapping.StructureMappingParameter) Parameter(cbit.vcell.model.Parameter) KineticsParameter(cbit.vcell.model.Kinetics.KineticsParameter) SpeciesContextSpecProxyParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecProxyParameter) ModelParameter(cbit.vcell.model.Model.ModelParameter) SpeciesContextSpecParameter(cbit.vcell.mapping.SpeciesContextSpec.SpeciesContextSpecParameter) ParticleProperties(cbit.vcell.math.ParticleProperties)

Example 28 with MembraneSubDomain

use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.

the class SimulationWarning method analyzeArea.

private static void analyzeArea(Simulation simulation, double timeStep, Map<MembraneSubDomain, List<DiffusionValue>> diffusionValuesMap, IssueContext issueContext, List<Issue> issueList) {
    double dx = simulation.getMeshSpecification().getDx(true);
    double dy = simulation.getMeshSpecification().getDy(true);
    double dz = simulation.getMeshSpecification().getDz(true);
    int dim = simulation.getMathDescription().getGeometry().getDimension();
    if (dim < 3) {
        throw new RuntimeException("suggested timestep analysis for smoldyn not implemented for " + dim + "D simulations");
    }
    double smallestQuadArea = Math.min(Math.min(dx * dy, dy * dz), dx * dz);
    // 1/2 for quad to triangle, 1/2 for smoother surface
    double averageTriangleSize = smallestQuadArea / 4;
    for (Entry<MembraneSubDomain, List<DiffusionValue>> entry : diffusionValuesMap.entrySet()) {
        MembraneSubDomain msd = entry.getKey();
        for (DiffusionValue dv : entry.getValue()) {
            double ad = averageTriangleSize / dv.value;
            boolean warn = (timeStep >= ad);
            if (lg.isDebugEnabled()) {
                lg.debug("average area " + averageTriangleSize + " diffusion " + dv.value + " time step limit " + ad + "timeStep " + timeStep + " -> warn " + warn);
            }
            if (warn) {
                DecimalFormat f = new DecimalFormat("0.##E0");
                String limit = f.format(ad);
                String m = "Solver time step " + timeStep + " may yield inaccurate results for given mesh size and diffusion coefficient of species " + dv.name + " on membrane " + msd.getName() + "; it is recommended that the time step not exceed " + limit;
                Issue i = new Issue(simulation, issueContext, IssueCategory.SMOLYDN_DIFFUSION, m, m, Severity.WARNING);
                issueList.add(i);
                if (lg.isDebugEnabled()) {
                    lg.debug("Warning " + m);
                }
            }
        }
    }
}
Also used : MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) Issue(org.vcell.util.Issue) DecimalFormat(java.text.DecimalFormat) ArrayList(java.util.ArrayList) List(java.util.List)

Example 29 with MembraneSubDomain

use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.

the class XmlReader method getMembraneSubDomain.

/**
 * This method returns a MembraneSubDomain object from a XML Element.
 * Creation date: (5/18/2001 4:23:30 PM)
 * @return cbit.vcell.math.MembraneSubDomain
 * @param param org.jdom.Element
 * @exception cbit.vcell.xml.XmlParseException The exception description.
 */
@SuppressWarnings("unchecked")
private MembraneSubDomain getMembraneSubDomain(Element param, MathDescription mathDesc) throws XmlParseException {
    // no need to do anything with the 'Name' attribute : constructor of MembraneSubDomain creates name from inside/outside compartmentSubDomains.
    // String msdName = unMangle( param.getAttributeValue(XMLTags.NameAttrTag) );
    // if ( msdName != null) {
    // }
    // get compartmentSubDomain references
    // inside
    String name = unMangle(param.getAttributeValue(XMLTags.InsideCompartmentTag));
    CompartmentSubDomain insideRef = (CompartmentSubDomain) mathDesc.getCompartmentSubDomain(name);
    if (insideRef == null) {
        throw new XmlParseException("The reference to the inside CompartmentSubDomain " + name + ", could not be resolved!");
    }
    // outside
    name = unMangle(param.getAttributeValue(XMLTags.OutsideCompartmentTag));
    CompartmentSubDomain outsideRef = (CompartmentSubDomain) mathDesc.getCompartmentSubDomain(name);
    if (outsideRef == null) {
        throw new XmlParseException("The reference to the outside CompartmentSubDomain " + name + ", could not be resolved!");
    }
    // *** create new Membrane SubDomain ***
    SubVolume insideSubVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(insideRef.getName());
    SubVolume outsideSubVolume = mathDesc.getGeometry().getGeometrySpec().getSubVolume(outsideRef.getName());
    SurfaceClass surfaceClass = mathDesc.getGeometry().getGeometrySurfaceDescription().getSurfaceClass(insideSubVolume, outsideSubVolume);
    MembraneSubDomain subDomain = new MembraneSubDomain(insideRef, outsideRef, surfaceClass.getName());
    transcribeComments(param, subDomain);
    // Process BoundaryConditions
    Iterator<Element> iterator = param.getChildren(XMLTags.BoundaryTypeTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        // create BoundaryConditionType
        String temp = tempelement.getAttributeValue(XMLTags.BoundaryTypeAttrTag);
        BoundaryConditionType bType = new BoundaryConditionType(temp);
        // Process Xm
        if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueXm)) {
            subDomain.setBoundaryConditionXm(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueXp)) {
            // Process Xp
            subDomain.setBoundaryConditionXp(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueYm)) {
            // Process Ym
            subDomain.setBoundaryConditionYm(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueYp)) {
            // Process Yp
            subDomain.setBoundaryConditionYp(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueZm)) {
            // Process Zm
            subDomain.setBoundaryConditionZm(bType);
        } else if (tempelement.getAttributeValue(XMLTags.BoundaryAttrTag).equalsIgnoreCase(XMLTags.BoundaryAttrValueZp)) {
            // Process Zp
            subDomain.setBoundaryConditionZp(bType);
        } else {
            // If not indentified throw an exception!!
            throw new XmlParseException("Unknown BoundaryConditionType: " + tempelement.getAttributeValue(XMLTags.BoundaryAttrTag));
        }
    }
    // Add OdeEquations
    iterator = param.getChildren(XMLTags.OdeEquationTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        OdeEquation odeEquation = getOdeEquation(tempElement, mathDesc);
        try {
            subDomain.addEquation(odeEquation);
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding an OdeEquation to a MembraneSubDomain!", e);
        }
    }
    // process PdeEquations
    iterator = param.getChildren(XMLTags.PdeEquationTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        try {
            subDomain.addEquation(getPdeEquation(tempElement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding an PdeEquation to the MembraneSubDomain " + name, e);
        }
    }
    // Add JumpConditions
    iterator = param.getChildren(XMLTags.JumpConditionTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempElement = (Element) iterator.next();
        try {
            subDomain.addJumpCondition(getJumpCondition(tempElement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a JumpCondition to a MembraneSubDomain!", e);
        }
    }
    // Add the FastSystem (if any)
    Element tempElement = param.getChild(XMLTags.FastSystemTag, vcNamespace);
    if (tempElement != null) {
        subDomain.setFastSystem(getFastSystem(tempElement, mathDesc));
    }
    // add MembraneRegionEquation
    iterator = param.getChildren(XMLTags.MembraneRegionEquationTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        tempElement = (Element) iterator.next();
        try {
            subDomain.addEquation(getMembraneRegionEquation(tempElement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a MembraneRegionEquation to a MEmbraneSubDomain!", e);
        }
    }
    iterator = param.getChildren(XMLTags.ParticleJumpProcessTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        try {
            subDomain.addParticleJumpProcess(getParticleJumpProcess(tempelement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a jump process to the MembraneSubDomain " + name, e);
        }
    }
    iterator = param.getChildren(XMLTags.ParticlePropertiesTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        try {
            subDomain.addParticleProperties(getParticleProperties(tempelement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding a jump process to the MembraneSubDomain " + name, e);
        }
    }
    // process ComputeNormal "equations"
    iterator = param.getChildren(XMLTags.ComputeNormalTag, vcNamespace).iterator();
    while (iterator.hasNext()) {
        Element tempelement = (Element) iterator.next();
        try {
            subDomain.addEquation(getComputeNormal(tempelement, mathDesc));
        } catch (MathException e) {
            e.printStackTrace();
            throw new XmlParseException("A MathException was fired when adding an ComputeNormal 'equation' to the MembraneSubDomain " + name, e);
        }
    }
    Element velElem = param.getChild(XMLTags.VelocityTag, vcNamespace);
    setMembraneSubdomainVelocity(velElem, XMLTags.XAttrTag, subDomain::setVelocityX);
    setMembraneSubdomainVelocity(velElem, XMLTags.YAttrTag, subDomain::setVelocityY);
    return subDomain;
}
Also used : MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) OdeEquation(cbit.vcell.math.OdeEquation) SurfaceClass(cbit.vcell.geometry.SurfaceClass) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubVolume(cbit.vcell.geometry.SubVolume) CompartmentSubVolume(cbit.vcell.geometry.CompartmentSubVolume) AnalyticSubVolume(cbit.vcell.geometry.AnalyticSubVolume) ImageSubVolume(cbit.vcell.geometry.ImageSubVolume) MathException(cbit.vcell.math.MathException) Element(org.jdom.Element) BoundaryConditionType(cbit.vcell.math.BoundaryConditionType)

Example 30 with MembraneSubDomain

use of cbit.vcell.math.MembraneSubDomain in project vcell by virtualcell.

the class SubdomainInfo method write.

public static void write(File file, MathDescription mathDesc) throws IOException, MathException {
    PrintWriter pw = null;
    try {
        pw = new PrintWriter(new FileWriter(file));
        pw.println("# " + VCML.CompartmentSubDomain + " name, handle");
        pw.println("# " + VCML.MembraneSubDomain + " name, inside compartment name, handle, outside compartment name, handle");
        Enumeration<SubDomain> subdomains = mathDesc.getSubDomains();
        while (subdomains.hasMoreElements()) {
            SubDomain sd = subdomains.nextElement();
            if (sd instanceof CompartmentSubDomain) {
                CompartmentSubDomain csd = (CompartmentSubDomain) sd;
                pw.println(VCML.CompartmentSubDomain + ", " + csd.getName() + ", " + mathDesc.getHandle(csd));
            } else if (sd instanceof MembraneSubDomain) {
                MembraneSubDomain msd = (MembraneSubDomain) sd;
                CompartmentSubDomain insideCompartment = msd.getInsideCompartment();
                CompartmentSubDomain outsideCompartment = msd.getOutsideCompartment();
                pw.println(VCML.MembraneSubDomain + ", " + msd.getName() + ", " + insideCompartment.getName() + ", " + mathDesc.getHandle(insideCompartment) + ", " + outsideCompartment.getName() + ", " + mathDesc.getHandle(outsideCompartment));
            }
        }
        pw.close();
    } finally {
        if (pw != null) {
            pw.close();
        }
    }
}
Also used : CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) SubDomain(cbit.vcell.math.SubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) MembraneSubDomain(cbit.vcell.math.MembraneSubDomain) CompartmentSubDomain(cbit.vcell.math.CompartmentSubDomain) FileWriter(java.io.FileWriter) PrintWriter(java.io.PrintWriter)

Aggregations

MembraneSubDomain (cbit.vcell.math.MembraneSubDomain)31 CompartmentSubDomain (cbit.vcell.math.CompartmentSubDomain)26 SubDomain (cbit.vcell.math.SubDomain)21 Expression (cbit.vcell.parser.Expression)13 ArrayList (java.util.ArrayList)13 MathDescription (cbit.vcell.math.MathDescription)12 SubVolume (cbit.vcell.geometry.SubVolume)10 ExpressionException (cbit.vcell.parser.ExpressionException)10 Variable (cbit.vcell.math.Variable)9 SurfaceClass (cbit.vcell.geometry.SurfaceClass)8 Constant (cbit.vcell.math.Constant)7 Equation (cbit.vcell.math.Equation)7 MathException (cbit.vcell.math.MathException)7 OdeEquation (cbit.vcell.math.OdeEquation)7 JumpCondition (cbit.vcell.math.JumpCondition)6 ParticleProperties (cbit.vcell.math.ParticleProperties)6 HashMap (java.util.HashMap)6 AnalyticSubVolume (cbit.vcell.geometry.AnalyticSubVolume)5 Function (cbit.vcell.math.Function)5 MemVariable (cbit.vcell.math.MemVariable)5