Skip to content

Commit a5d6100

Browse files
tomolearytomoleary
authored
Initial code (#1)
* testing readme * testing readme * testing readme * testing readme * adding surrogate construction files * updating
1 parent a068327 commit a5d6100

31 files changed

+4565
-3
lines changed

‎.gitignore

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
*.DS_Store
2+
13
# Byte-compiled / optimized / DLL files
24
__pycache__/
35
*.py[cod]

‎README.md

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
1-
Derivative Informed Neural Operator
2-
An Efficient Framework for High-Dimensional Parametric Derivative Learning
1+
Derivative Informed Neural Operator
32

43
_____ ___ ___
54
/ /::\ ___ /__/\ / /\
@@ -13,4 +12,4 @@
1312
\__\/ \__\/ \ \:\ \ \::/
1413
\__\/ \__\/
1514

16-
15+
An Efficient Framework for High-Dimensional Parametric Derivative Learning

‎applications/confusion/confusionLinearObservable.py

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
# This file is part of the dino package
2+
#
3+
# dino is free software: you can redistribute it and/or modify
4+
# it under the terms of the GNU Lesser General Public License as published by
5+
# the Free Software Foundation, either version 2 of the License, or any later version.
6+
#
7+
# dino is distributed in the hope that it will be useful,
8+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
9+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10+
# GNU Lesser General Public License for more details.
11+
#
12+
# You should have received a copy of the GNU Lesser General Public License
13+
# If not, see <http://www.gnu.org/licenses/>.
14+
#
15+
# Author: Tom O'Leary-Roseberry
16+
# Contact: [email protected]
17+
18+
import dolfin as dl
19+
import numpy as np
20+
import ufl
21+
22+
path_to_hippylib = '../../hippylib/'
23+
import sys, os
24+
sys.path.append( os.environ.get('HIPPYLIB_PATH',path_to_hippylib))
25+
from hippylib import *
26+
27+
sys.path.append( os.environ.get('HIPPYFLOW_PATH'))
28+
from hippyflow import LinearStateObservable
29+
30+
def confusion_linear_observable(mesh,nx_targets= 10, ny_targets = 5,use_krylov = False,output_folder ='poisson_setup/',\
31+
verbose = False,seed = 0):
32+
class confusion_varf:
33+
def __init__(self,Vh,save_fields = False,output_folder = 'confusion_setup/'):
34+
'''
35+
36+
'''
37+
self.Vh = Vh
38+
# Gaussian blob for the right hand side
39+
self.f = dl.interpolate( dl.Expression('max(0.5,exp(-25*(pow(x[0]-0.7,2) + pow(x[1]-0.7,2))))',degree=5), Vh[STATE])
40+
# Compute velocity field
41+
self.v = self.computeVelocityField(Vh[STATE].mesh())
42+
if save_fields:
43+
if not os.path.isdir(output_folder):
44+
os.mkdir(output_folder)
45+
f_pvd = dl.File(output_folder+'f_blob.pvd')
46+
f_pvd << self.f
47+
v_pvd = dl.File(output_folder+'v_sol.pvd')
48+
v_pvd << self.v
49+
# Constant coefficients for the PDE
50+
self.c = dl.Constant(1.0)
51+
self.k = dl.Constant(0.01)
52+
53+
54+
def computeVelocityField(self,mesh):
55+
def v_boundary(x,on_boundary):
56+
return on_boundary
57+
58+
def q_boundary(x,on_boundary):
59+
return x[0] < dl.DOLFIN_EPS and x[1] < dl.DOLFIN_EPS
60+
61+
Xh = dl.VectorFunctionSpace(mesh,'Lagrange', 2)
62+
Wh = dl.FunctionSpace(mesh, 'Lagrange', 1)
63+
mixed_element = ufl.MixedElement([Xh.ufl_element(), Wh.ufl_element()])
64+
XW = dl.FunctionSpace(mesh, mixed_element)
65+
66+
Re = dl.Constant(1e2)
67+
68+
g = dl.Expression(('0.0','(x[0] < 1e-14) - (x[0] > 1 - 1e-14)'), degree=1)
69+
bc1 = dl.DirichletBC(XW.sub(0), g, v_boundary)
70+
bc2 = dl.DirichletBC(XW.sub(1), dl.Constant(0), q_boundary, 'pointwise')
71+
bcs = [bc1, bc2]
72+
73+
vq = dl.Function(XW)
74+
(v,q) = ufl.split(vq)
75+
(v_test, q_test) = dl.TestFunctions (XW)
76+
77+
def strain(v):
78+
return ufl.sym(ufl.grad(v))
79+
80+
F = ( (2./Re)*ufl.inner(strain(v),strain(v_test))+ ufl.inner (ufl.nabla_grad(v)*v, v_test)
81+
- (q * ufl.div(v_test)) + ( ufl.div(v) * q_test) ) * ufl.dx
82+
83+
dl.solve(F == 0, vq, bcs, solver_parameters={"newton_solver":
84+
{"relative_tolerance":1e-4, "maximum_iterations":150}})
85+
return vq.split()[0]
86+
87+
88+
def __call__(self,u,m,p):
89+
'''
90+
Return the variational form of the PDE
91+
Inputs
92+
-u: state variable
93+
-m: model parameter
94+
-p: adjoint variable
95+
Outputs:
96+
Variational form of the PDE
97+
'''
98+
h = dl.CellDiameter(Vh[STATE].mesh())
99+
v_norm = dl.sqrt( dl.dot(self.v,self.v) + 1e-6)
100+
101+
return (h/v_norm)*dl.dot( self.v, dl.nabla_grad(u)) * dl.dot( self.v, dl.nabla_grad(p)) * dl.dx\
102+
+ self.k*dl.inner(dl.nabla_grad(u), dl.nabla_grad(p))*dl.dx \
103+
+ dl.inner(dl.nabla_grad(u), self.v*p)*dl.dx \
104+
+ self.c*dl.exp(m)*u*u*u*p*dl.dx \
105+
- self.f*p*dl.dx
106+
107+
108+
109+
def u_boundary(x, on_boundary):
110+
return on_boundary
111+
112+
########################################################################
113+
114+
# Construct the linear observable
115+
116+
# Define the function spaces
117+
Vh1 = dl.FunctionSpace(mesh, 'Lagrange', 1)
118+
Vh = [Vh1, Vh1, Vh1]
119+
if verbose:
120+
print( "Number of dofs: STATE={0}, PARAMETER={1}, ADJOINT={2}".format(Vh[STATE].dim(), Vh[PARAMETER].dim(), Vh[ADJOINT].dim()) )
121+
122+
########################################################################
123+
# Construct the linear observable
124+
# Targets only on the bottom
125+
#Targets only on the bottom
126+
x_targets = np.linspace(0.1,0.9,nx_targets)
127+
y_targets = np.linspace(0.1,0.5,ny_targets)
128+
targets = []
129+
for xi in x_targets:
130+
for yi in y_targets:
131+
targets.append((xi,yi))
132+
targets = np.array(targets)
133+
ntargets = targets.shape[0]
134+
135+
if verbose:
136+
print( "Number of observation points: {0}".format(targets.shape[0]) )
137+
138+
# Define Dirichp.et boundary conditions
139+
u_bdr = dl.Constant(0.0)
140+
u_bdr0 = dl.Constant(0.0)
141+
bc = dl.DirichletBC(Vh[STATE], u_bdr, u_boundary)
142+
bc0 = dl.DirichletBC(Vh[STATE], u_bdr0, u_boundary)
143+
144+
145+
m0 = dl.interpolate(dl.Constant(0.0), Vh[PARAMETER]).vector()
146+
param_dimension = m0.get_local().shape[0]
147+
m0.set_local(np.random.randn(param_dimension))
148+
149+
varf_handler = confusion_varf(Vh, output_folder = output_folder)
150+
151+
pde = PDEVariationalProblem(Vh, varf_handler, bc, bc0, is_fwd_linear=False)
152+
if use_krylov:
153+
pde.solver = PETScKrylovSolver(mesh.mpi_comm(), "gmres", amg_method())
154+
pde.solver_fwd_inc = PETScKrylovSolver(mesh.mpi_comm(), "gmres", amg_method())
155+
pde.solver_adj_inc = PETScKrylovSolver(mesh.mpi_comm(), "gmres", amg_method())
156+
pde.solver.parameters["relative_tolerance"] = 1e-6
157+
pde.solver.parameters["absolute_tolerance"] = 1e-8
158+
pde.solver.parameters["maximum_iterations"] = 500
159+
pde.solver_fwd_inc.parameters = pde.solver.parameters
160+
pde.solver_fwd_inc.parameters["relative_tolerance"] = 1e-2
161+
pde.solver_adj_inc.parameters = pde.solver.parameters
162+
pde.solver_adj_inc.parameters["relative_tolerance"] = 1e-2
163+
164+
B = assemblePointwiseObservation(Vh[STATE], targets)
165+
166+
observable = LinearStateObservable(pde,B)
167+
168+
observable._targets_ = targets
169+
170+
return observable
171+
172+

‎applications/confusion/confusionModelSettings.py

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
# This file is part of the dino package
2+
#
3+
# dino is free software: you can redistribute it and/or modify
4+
# it under the terms of the GNU Lesser General Public License as published by
5+
# the Free Software Foundation, either version 2 of the License, or any later version.
6+
#
7+
# dino is distributed in the hope that it will be useful,
8+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
9+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10+
# GNU Lesser General Public License for more details.
11+
#
12+
# You should have received a copy of the GNU Lesser General Public License
13+
# If not, see <http://www.gnu.org/licenses/>.
14+
#
15+
# Author: Tom O'Leary-Roseberry
16+
# Contact: [email protected]
17+
18+
import math
19+
20+
def confusion_problem_settings(settings = {}):
21+
"""
22+
"""
23+
24+
settings['ntargets'] = 50
25+
settings['nx_targets'] = 10
26+
settings['ny_targets'] = 5
27+
settings['gamma'] = 0.1
28+
settings['delta'] = 1.0
29+
settings['nx'] = 64
30+
settings['ny'] = 64
31+
settings['ndim'] = 2
32+
settings['jacobian_full_rank'] = 50
33+
settings['formulation'] = 'confusion'
34+
settings['rel_noise'] = 0.01
35+
36+
# For prior anisotropic tensor
37+
settings['prior_theta0'] = 2.0
38+
settings['prior_theta1'] = 0.5
39+
settings['prior_alpha'] = math.pi/4
40+
41+
# For sampling
42+
settings['seed'] = 1
43+
44+
45+
return settings

‎applications/confusion/confusionModelUtilities.py

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
# This file is part of the dino package
2+
#
3+
# dino is free software: you can redistribute it and/or modify
4+
# it under the terms of the GNU Lesser General Public License as published by
5+
# the Free Software Foundation, either version 2 of the License, or any later version.
6+
#
7+
# dino is distributed in the hope that it will be useful,
8+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
9+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10+
# GNU Lesser General Public License for more details.
11+
#
12+
# You should have received a copy of the GNU Lesser General Public License
13+
# If not, see <http://www.gnu.org/licenses/>.
14+
#
15+
# Author: Tom O'Leary-Roseberry
16+
# Contact: [email protected]
17+
18+
import dolfin as dl
19+
import ufl
20+
import math
21+
import numpy as np
22+
import matplotlib.pyplot as plt
23+
import sys
24+
import os
25+
sys.path.append( os.environ.get('HIPPYLIB_PATH', "../") )
26+
import hippylib as hp
27+
28+
sys.path.append( os.environ.get('HIPPYFLOW_PATH'))
29+
import hippyflow as hf
30+
31+
import logging
32+
logging.getLogger('FFC').setLevel(logging.WARNING)
33+
logging.getLogger('UFL').setLevel(logging.WARNING)
34+
dl.set_log_active(False)
35+
36+
np.random.seed(seed=0)
37+
38+
from confusionModelSettings import confusion_problem_settings
39+
from confusionLinearObservable import confusion_linear_observable
40+
41+
42+
def confusion_model_wrapper(settings = confusion_problem_settings()):
43+
# Set up the mesh, finite element spaces, and PDE forward problem
44+
45+
ndim = settings['ndim']
46+
nx = settings['nx']
47+
ny = settings['ny']
48+
mesh = dl.UnitSquareMesh(nx, ny)
49+
50+
nx_targets = settings['nx_targets']
51+
ny_targets = settings['ny_targets']
52+
53+
observable_kwargs = {'nx_targets':nx_targets,'ny_targets':ny_targets}
54+
observable = confusion_linear_observable(mesh,**observable_kwargs)
55+
56+
pde = observable.problem
57+
Vh = pde.Vh
58+
# Set up the prior
59+
60+
gamma = settings['gamma']
61+
delta = settings['delta']
62+
63+
theta0 = settings['prior_theta0']
64+
theta1 = settings['prior_theta1']
65+
alpha = settings['prior_alpha']
66+
67+
anis_diff = dl.CompiledExpression(hp.ExpressionModule.AnisTensor2D(), degree = 1)
68+
anis_diff.set(theta0, theta1, alpha)
69+
70+
# Set up the prior
71+
prior = hp.BiLaplacianPrior(Vh[hp.PARAMETER], gamma, delta, anis_diff, robin_bc=True)
72+
print("Prior regularization: (delta_x - gamma*Laplacian)^order: delta={0}, gamma={1}, order={2}".format(delta, gamma,2))
73+
74+
targets = observable._targets_
75+
76+
ntargets = len(targets)
77+
78+
print( "Number of observation points: {0}".format(ntargets) )
79+
80+
misfit = hp.PointwiseStateObservation(Vh[hp.STATE], targets)
81+
82+
model = hp.Model(pde, prior, misfit)
83+
84+
mw_settings = hf.hippylibModelWrapperSettings()
85+
mw_settings['rel_noise'] = settings['rel_noise']
86+
mw_settings['seed'] = settings['seed']
87+
88+
modelwrapper = hf.hippylibModelWrapper(model,settings = mw_settings)
89+
modelwrapper.targets = targets
90+
# Setup the inverse problem
91+
modelwrapper.setUpInverseProblem()
92+
93+
return modelwrapper
94+

0 commit comments

Comments
 (0)