Skip to content

Commit

Permalink
Enabled using multiple regModels (#629)
Browse files Browse the repository at this point in the history 
* Enabled using multiple regModels.

* Fixed a warning for location obj.

* Fixed a regModel bug for forward.
  • Loading branch information
@friedenhe
friedenhe committed Apr 14, 2024
1 parent 528b283 commit be232cb
Show file tree
Hide file tree
Showing 20 changed files with 558 additions and 450 deletions.
22 changes: 18 additions & 4 deletions dafoam/mphys/mphys_dafoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -664,7 +664,8 @@ def setup(self):
distributed = designVariables[dvName]["distributed"]
self.add_input(dvName, distributed=distributed, shape_by_conn=True, tags=["mphys_coupling"])
elif dvType == "RegPar":
nParameters = self.DASolver.solver.getNRegressionParameters()
modelName = designVariables[dvName]["modelName"]
nParameters = self.DASolver.solver.getNRegressionParameters(modelName.encode())
self.add_input(dvName, distributed=False, shape=nParameters, tags=["mphys_coupling"])
else:
raise AnalysisError("designVarType %s not supported! " % dvType)
Expand Down Expand Up @@ -915,7 +916,11 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
product = np.zeros_like(parameters)
seeds = d_residuals["%s_states" % self.discipline]

DASolver.solverAD.calcdRdRegParTPsiAD(volCoords, states, parameters, seeds, product)
modelName = designVariables[inputName]["modelName"]

DASolver.solverAD.calcdRdRegParTPsiAD(
volCoords, states, parameters, seeds, modelName.encode(), product
)

# reduce to make sure all procs get the same product
product = self.comm.allreduce(product, op=MPI.SUM)
Expand Down Expand Up @@ -1238,7 +1243,8 @@ def setup(self):
distributed = designVariables[dvName]["distributed"]
self.add_input(dvName, distributed=distributed, shape_by_conn=True, tags=["mphys_coupling"])
elif dvType == "RegPar":
nParameters = self.DASolver.solver.getNRegressionParameters()
modelName = designVariables[dvName]["modelName"]
nParameters = self.DASolver.solver.getNRegressionParameters(modelName.encode())
self.add_input(dvName, distributed=False, shape=nParameters, tags=["mphys_coupling"])
else:
raise AnalysisError("designVarType %s not supported! " % dvType)
Expand Down Expand Up @@ -1456,8 +1462,16 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
parameters = inputs[inputName]
product = np.zeros_like(parameters)

modelName = designVariables[inputName]["modelName"]

DASolver.solverAD.calcdFdRegParAD(
volCoords, states, parameters, objFuncName.encode(), inputName.encode(), product
volCoords,
states,
parameters,
objFuncName.encode(),
inputName.encode(),
modelName.encode(),
product,
)

d_inputs[inputName] += product * fBar
Expand Down
80 changes: 54 additions & 26 deletions dafoam/pyDAFoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -514,21 +514,43 @@ def __init__(self):
self.useConstrainHbyA = False

## parameters for regression models
## we support defining multiple regression models. Each regression model can have only one output
## but it can have multiple input features. Refer to src/adjoint/DARegression/DARegression.C for
## a full list of supported input features. There are two supported regression model types:
## neural network and radial basis function. We can shift and scale the inputs and outputs
## we can also prescribe a default output value. The default output will be used in resetting
## when they are nan, inf, or out of the prescribe upper and lower bounds
self.regressionModel = {
"active": False,
"modelType": "neuralNetwork",
"inputNames": ["None"],
"outputName": "None",
"hiddenLayerNeurons": [0],
"inputShift": [0.0],
"inputScale": [1.0],
"outputShift": 0.0,
"outputScale": 1.0,
"outputUpperBound": 1e8,
"outputLowerBound": -1e8,
"activationFunction": "sigmoid",
"printInputInfo": True,
"defaultOutputValue": 0.0,
# "model1": {
# "modelType": "neuralNetwork",
# "inputNames": ["PoD", "VoS", "PSoSS", "chiSA"],
# "outputName": "betaFINuTilda",
# "hiddenLayerNeurons": [20, 20],
# "inputShift": [0.0],
# "inputScale": [1.0],
# "outputShift": 1.0,
# "outputScale": 1.0,
# "outputUpperBound": 1e1,
# "outputLowerBound": -1e1,
# "activationFunction": "sigmoid", # other options are relu and tanh
# "printInputInfo": True,
# "defaultOutputValue": 1.0,
# },
# "model2": {
# "modelType": "radialBasisFunction",
# "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
# "outputName": "betaFIOmega",
# "nRBFs": 50,
# "inputShift": [0.0],
# "inputScale": [1.0],
# "outputShift": 1.0,
# "outputScale": 1.0,
# "outputUpperBound": 1e1,
# "outputLowerBound": -1e1,
# "printInputInfo": True,
# "defaultOutputValue": 1.0,
# },
}

## whether to use averaged states
Expand All @@ -539,10 +561,7 @@ def __init__(self):
## we use no time steps for averaging. 0.8 means we use the last 20% of the time step for averaging.
## We usually don't use 0 because the flow will need some spin up time, so using the spin-up
## flow field for meanStates will be slightly inaccurate.
self.useMeanStates = {
"active": False,
"start": 0.5
}
self.useMeanStates = {"active": False, "start": 0.5}

# *********************************************************************************************
# ************************************ Advance Options ****************************************
Expand Down Expand Up @@ -2181,7 +2200,7 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin
else:
self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivSeq[i]

def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accumulateTotal=False):
def calcTotalDerivsRegPar(self, objFuncName, designVarName, modelName, dFScaling=1.0, accumulateTotal=False):

nDVs = 0
parameters = None
Expand All @@ -2195,7 +2214,7 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accum
nDVs = len(self.internalDV[designVarName]["init"])
parameters = self.internalDV[designVarName]["value"]

nParameters = self.solver.getNRegressionParameters()
nParameters = self.getNRegressionParameters(modelName)
if nDVs != nParameters:
raise Error("number of parameters not valid!")

Expand All @@ -2207,12 +2226,12 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accum

# calc dFdRegPar
self.solverAD.calcdFdRegParAD(
xvArray, wArray, parameters, objFuncName.encode(), designVarName.encode(), dFdRegPar
xvArray, wArray, parameters, objFuncName.encode(), designVarName.encode(), modelName.encode(), dFdRegPar
)
dFdRegPar *= dFScaling

# calculate dRdFieldT*Psi and save it to totalDeriv
self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, totalDerivArray)
self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, modelName.encode(), totalDerivArray)
# all reduce because parameters is a global DV
totalDerivArray = self.comm.allreduce(totalDerivArray, op=MPI.SUM)

Expand Down Expand Up @@ -2437,7 +2456,8 @@ def solveAdjointUnsteady(self):
fieldType = designVarDict[designVarName]["fieldType"]
self.calcTotalDerivsField(objFuncName, designVarName, fieldType, dFScaling, True)
elif designVarDict[designVarName]["designVarType"] == "RegPar":
self.calcTotalDerivsRegPar(objFuncName, designVarName, dFScaling, True)
modelName = designVarDict[designVarName]["modelName"]
self.calcTotalDerivsRegPar(objFuncName, designVarName, modelName, dFScaling, True)
else:
raise Error("designVarType not valid!")

Expand Down Expand Up @@ -2736,7 +2756,8 @@ def solveAdjoint(self):
# loop over all objectives
for objFuncName in objFuncDict:
if objFuncName in self.objFuncNames4Adj:
self.calcTotalDerivsRegPar(objFuncName, designVarName)
modelName = designVarDict[designVarName]["modelName"]
self.calcTotalDerivsRegPar(objFuncName, designVarName, modelName)
else:
raise Error("For Field design variable type, we only support useAD->mode=reverse")
else:
Expand Down Expand Up @@ -4013,7 +4034,7 @@ def setThermal(self, thermalArray):
if self.getOption("useAD")["mode"] in ["forward", "reverse"]:
self.solverAD.setThermal(thermalArray)

def setRegressionParameter(self, idx, val):
def setRegressionParameter(self, modelName, idx, val):
"""
Update the regression parameters
Expand All @@ -4026,9 +4047,16 @@ def setRegressionParameter(self, idx, val):
the parameter value to set
"""

self.solver.setRegressionParameter(idx, val)
self.solver.setRegressionParameter(modelName.encode(), idx, val)
if self.getOption("useAD")["mode"] in ["forward", "reverse"]:
self.solverAD.setRegressionParameter(idx, val)
self.solverAD.setRegressionParameter(modelName.encode(), idx, val)

def getNRegressionParameters(self, modelName):
"""
Get the number of regression model parameters
"""
nParameters = self.solver.getNRegressionParameters(modelName.encode())
return nParameters

def setFieldValue4GlobalCellI(self, fieldName, val, globalCellI, compI=0):
"""
Expand Down
2 changes: 1 addition & 1 deletion src/adjoint/DAObjFunc/DAObjFuncLocation.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ DAObjFuncLocation::DAObjFuncLocation(
// we assume the patchI and faceI do not change during the optimization
// otherwise, we should use maxRadiusKS instead
scalar maxR = -100000;
label maxRPatchI, maxRFaceI;
label maxRPatchI = -999, maxRFaceI = -999;
forAll(objFuncFaceSources_, idxI)
{
const label& objFuncFaceI = objFuncFaceSources_[idxI];
Expand Down

0 comments on commit be232cb

Please sign in to comment.