I gave the Migration Guide a read.  Thanks for sending me that!  I do suspect once I get to the bottom of this something from that guide will be what I did wrong since my custom functions were using .data in 0.3 and 1.3 versions.  However, when making the minimal example my custom functions were not needed to replicate this problem.  The following code also shows that the hand calculated grad does not match up with the autograd calculated one in 1.3.1.  I made similar code in 0.3.0 and it also does not match there.
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F

import torchvision.transforms as transforms
from torchvision import datasets

class Small_Net(nn.Module):
    def __init__(self, device):
        super(Small_Net, self).__init__()
        self.layer1 = nn.Conv2d(1, 40, kernel_size=28, stride=1, padding=0, bias=True)
        self.layer2 = nn.Conv2d(40, 10, kernel_size=1, stride=1, padding=0, bias=True)
        self.layer1.to(device)
        self.layer2.to(device)

    def forward(self, x):        
        out2 = self.layer1(x).sigmoid()
        self.midOut = out2.detach().clone()
        out3 = self.layer2(out2)            
        out6 = out3.view(out3.size(0), -1)
        return out6

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
net = Small_Net(device)  
net.to(device)
net = torch.nn.DataParallel(net, device_ids=range(torch.cuda.device_count())).to(device)
start_epoch = 0
trainloader = torch.utils.data.DataLoader(datasets.MNIST('../mnist_data', 
                                                          download=True, 
                                                          train=True,
                                                          transform=transforms.Compose([
                                                              transforms.ToTensor(), # first, convert image to PyTorch tensor
                                                              transforms.Normalize((0.1307,), (0.3081,)) # normalize inputs
                                                          ])), 
                                           batch_size=128, 
                                           shuffle=True)
criterion = nn.CrossEntropyLoss()             
def train(epoch, mode):
    net.train()
    optimizer = optim.SGD(filter(lambda p: p.requires_grad, net.parameters()), lr=0.1, momentum=0, weight_decay=0)  # 5e-4)
    for batch_idx, (inputs, targets) in enumerate(trainloader):
        batchSize, outputMaps, x, y = inputs.size()
        inputs, targets = inputs.to(device), targets.to(device)  # GPU settings        
        optimizer.zero_grad()
        output = net(inputs)
        loss = criterion(output, targets)
        loss.backward()    
        X = output
        y = targets
        m = y.shape[0]
        grad = F.softmax(X,1)
        grad[range(m),y] -= 1
        grad = grad/m
        print('is same?')        
        print(str((net.module.layer1(inputs).sigmoid()[:,0,0,0] * grad[:,0]).sum().tolist()))
        print((net.module.midOut[:,0,0,0] * grad[:,0]).sum().tolist())
        print(str((net.module.layer2.weight.grad[0][0][0][0]).tolist()))
        del inputs, targets, output
        optimizer.step() 
    
train(0, 'n')

Example output:

is same?
-0.0011598279234021902
-0.0011598279234021902
-0.0011598295532166958

Edit: with double precision it does agree much more but still not the same.  Is this just the best that can be expected because of what autograd is doing under the hood?

is same?
-0.00020605623917282047
-0.00020605623917282047
-0.00020605623917282427